diff --git a/examples/WIMP.rml b/examples/WIMP.rml index c944a6a..abdf50e 100644 --- a/examples/WIMP.rml +++ b/examples/WIMP.rml @@ -10,7 +10,7 @@ - + diff --git a/examples/WIMP_compound_1.rml b/examples/WIMP_compound_1.rml new file mode 100644 index 0000000..6762a4e --- /dev/null +++ b/examples/WIMP_compound_1.rml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/examples/WIMP_compound_2.rml b/examples/WIMP_compound_2.rml new file mode 100644 index 0000000..4501b89 --- /dev/null +++ b/examples/WIMP_compound_2.rml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/inc/TRestWimpNucleus.h b/inc/TRestWimpNucleus.h index de2267a..23f8a83 100644 --- a/inc/TRestWimpNucleus.h +++ b/inc/TRestWimpNucleus.h @@ -38,8 +38,11 @@ class TRestWimpNucleus { Int_t fZnum; /// Abundance, in mass percentage Double_t fAbundance; + /// Abundance, in mole (or volume) + Double_t fAbundanceMol; void PrintNucleus(); + int GetStechiometricFactorFromCompound(const std::string& compound); // Constructor TRestWimpNucleus(); @@ -47,7 +50,7 @@ class TRestWimpNucleus { // Destructor virtual ~TRestWimpNucleus(); - ClassDef(TRestWimpNucleus, 1); + ClassDef(TRestWimpNucleus, 2); }; #endif diff --git a/inc/TRestWimpUtils.h b/inc/TRestWimpUtils.h index 38f8030..e26c46c 100644 --- a/inc/TRestWimpUtils.h +++ b/inc/TRestWimpUtils.h @@ -21,6 +21,7 @@ *************************************************************************/ #include +#include #ifndef RestCore_TRestWimpUtils #define RestCore_TRestWimpUtils @@ -54,7 +55,7 @@ const double GetRecoilRate(const double wimpMass, const double crossSection, con const double Anum, const double vLab, const double vRMS, const double vEscape, const double wimpDensity, const double abundance); const double GetQuenchingFactor(const double recoilEnergy, const double Anum, const double Znum); - +std::map ParseChemicalCompound(const std::string& compound); } // namespace TRestWimpUtils #endif diff --git a/src/TRestWimpNucleus.cxx b/src/TRestWimpNucleus.cxx index 48b1fd7..567c519 100644 --- a/src/TRestWimpNucleus.cxx +++ b/src/TRestWimpNucleus.cxx @@ -42,10 +42,17 @@ #include "TRestWimpNucleus.h" #include "TRestMetadata.h" +#include "TRestWimpUtils.h" ClassImp(TRestWimpNucleus); -TRestWimpNucleus::TRestWimpNucleus() {} +TRestWimpNucleus::TRestWimpNucleus() { + fNucleusName = ""; + fAnum = 0; + fZnum = 0; + fAbundance = 0; + fAbundanceMol = 0; +} TRestWimpNucleus::~TRestWimpNucleus() {} @@ -55,5 +62,28 @@ void TRestWimpNucleus::PrintNucleus() { RESTMetadata << "Atomic number " << fAnum << RESTendl; RESTMetadata << "Number of protons " << fZnum << RESTendl; RESTMetadata << "Abundance " << fAbundance << RESTendl; + RESTMetadata << "Abundance (mol) " << fAbundanceMol << RESTendl; RESTMetadata << "-----------------------------" << RESTendl; } + +/////////////////////////////////////////////// +/// \brief Get the stechiometric factor of this nucleus in a given compound. +/// +/// \param compound The compound to be parsed. See +/// TRestWimpUtils::ParseChemicalCompound for the compound format. +/// +int TRestWimpNucleus::GetStechiometricFactorFromCompound(const std::string& compound) { + auto elementMap = TRestWimpUtils::ParseChemicalCompound(compound); + + int stechiometricFactor = 0; + for (auto& pair : elementMap) { + if (pair.first == fNucleusName.Data()) { + stechiometricFactor = pair.second; + break; + } + } + if (stechiometricFactor == 0) + RESTWarning << "No nucleus " << fNucleusName.Data() << " founnd in compound " << compound << RESTendl; + + return stechiometricFactor; +} diff --git a/src/TRestWimpSensitivity.cxx b/src/TRestWimpSensitivity.cxx index 7987e6c..acc2dfd 100644 --- a/src/TRestWimpSensitivity.cxx +++ b/src/TRestWimpSensitivity.cxx @@ -37,6 +37,8 @@ /// - *anum* : Atomic number of the nucleus /// - *znum* : Number of protons in the nucleus /// - *abundance* : Mass percentage of the nucleus in the target material. +/// - *abundanceInMol* : Mol (or volume) percentage of the nucleus in the +/// target material. /// - *wimpDensity* : WIMP density in the DM halo in GeV/cm3 /// - *labVelocity* : WIMP velocity in the lab (Earth) frame reference /// in km/s @@ -68,6 +70,12 @@ /// /// /// \endcode +/// Other ways to define the target material is through the use of +/// compounds. Also, the abundances can be given in mol (or volume) using +/// the parameter *abundanceInMol* instead of *abundance*. Examples for an +/// Ar+Isobutane mixture at 99% in volume can be found inside the files +/// 'REST_PATH/source/libraries/wimp/examples/WIMP_compound_1.rml' and +/// 'REST_PATH/source/libraries/wimp/examples/WIMP_compound_2.rml'. /// /// Besides target material elements, the other parameters are optional, /// and if not provided they will take their default values. @@ -161,23 +169,123 @@ void TRestWimpSensitivity::InitFromConfigFile() { /// void TRestWimpSensitivity::ReadNuclei() { fNuclei.clear(); + bool anyAbundanceGivenInMol = false; + // Read nuclei (standalone given) elements TiXmlElement* ElementDef = GetElement("addElement"); while (ElementDef) { TRestWimpNucleus nucleus; nucleus.fNucleusName = GetFieldValue("nucleusName", ElementDef); nucleus.fAnum = StringToDouble(GetFieldValue("anum", ElementDef)); nucleus.fZnum = StringToInteger(GetFieldValue("znum", ElementDef)); - nucleus.fAbundance = StringToDouble(GetFieldValue("abundance", ElementDef)); + + std::string el = + !ElementDef->Attribute("abundance") ? "Not defined" : ElementDef->Attribute("abundance"); + if (!(el.empty() || el == "Not defined")) nucleus.fAbundance = StringToDouble(el); + + el = !ElementDef->Attribute("abundanceInMol") ? "Not defined" + : ElementDef->Attribute("abundanceInMol"); + if (!(el.empty() || el == "Not defined")) { + nucleus.fAbundanceMol = StringToDouble(el); + anyAbundanceGivenInMol = true; + } + + if (nucleus.fAbundance == 0 || nucleus.fAbundanceMol == 0) { + if (nucleus.fAbundance == 0) + nucleus.fAbundance = nucleus.fAbundanceMol * nucleus.fAnum; + else if (nucleus.fAbundanceMol == 0) + nucleus.fAbundanceMol = nucleus.fAbundance / nucleus.fAnum; + else + RESTError << "abundance or abundanceInMol not defined for nucleus " << nucleus.fNucleusName + << RESTendl; + } fNuclei.emplace_back(nucleus); ElementDef = GetNextElement(ElementDef); } + + // Read nuclei (compound form given) elements + TiXmlElement* CompoundDef = GetElement("addCompound"); + while (CompoundDef) { + bool compoundAbundanceGivenInMol = false; + std::string compoundName = GetFieldValue("compoundName", CompoundDef); + double compoundAbundance = 0, compoundAbundanceInMol = 0; + double totalMolMass = 0; + + std::string el = + !CompoundDef->Attribute("abundance") ? "Not defined" : CompoundDef->Attribute("abundance"); + if (!(el.empty() || el == "Not defined")) compoundAbundance = StringToDouble(el); + + el = !CompoundDef->Attribute("abundanceInMol") ? "Not defined" + : CompoundDef->Attribute("abundanceInMol"); + if (!(el.empty() || el == "Not defined")) { + compoundAbundanceInMol = StringToDouble(el); + compoundAbundanceGivenInMol = true; + anyAbundanceGivenInMol = true; + } + + if (compoundAbundance == 0 && compoundAbundanceInMol == 0) { + RESTWarning << "abundance or abundanceInMol not defined for compound " << compoundName + << ". Setting its abundanceInMol to 1 " << RESTendl; + compoundAbundanceInMol = 1; + } + + // Read nuclei (inside compound) elements + TiXmlElement* ElementDef = GetElement("addElement", CompoundDef); + int i = 0; + while (ElementDef) { + i++; + TRestWimpNucleus nucleus; + nucleus.fNucleusName = GetFieldValue("nucleusName", ElementDef); + nucleus.fAnum = StringToDouble(GetFieldValue("anum", ElementDef)); + nucleus.fZnum = StringToInteger(GetFieldValue("znum", ElementDef)); + totalMolMass += nucleus.fAnum * nucleus.GetStechiometricFactorFromCompound(compoundName); + + fNuclei.emplace_back(nucleus); + ElementDef = GetNextElement(ElementDef); + } + if (compoundAbundanceGivenInMol) + compoundAbundance = compoundAbundanceInMol * totalMolMass; + else + compoundAbundanceInMol = compoundAbundance / totalMolMass; + // Set the compound abundance to all nuclei elements inside the compound + for (auto it = fNuclei.end() - i; it != fNuclei.end(); it++) { + auto& nucleus = *it; + int stechiometricFactor = nucleus.GetStechiometricFactorFromCompound(compoundName); + nucleus.fAbundanceMol = compoundAbundanceInMol * stechiometricFactor; + nucleus.fAbundance = nucleus.fAbundanceMol * nucleus.fAnum; + } + + CompoundDef = GetNextElement(CompoundDef); + } + + // Merge the repeated nuclei (same name, anum and znum) by summing their abundances + std::map, TRestWimpNucleus> uniqueNucleiMap; + for (const auto& nucleus : fNuclei) { + auto key = std::make_tuple(nucleus.fNucleusName, nucleus.fAnum, nucleus.fZnum); + if (uniqueNucleiMap.find(key) != uniqueNucleiMap.end()) { + uniqueNucleiMap[key].fAbundance += nucleus.fAbundance; + uniqueNucleiMap[key].fAbundanceMol += nucleus.fAbundanceMol; + } else + uniqueNucleiMap[key] = nucleus; + } + fNuclei.clear(); + for (const auto& entry : uniqueNucleiMap) fNuclei.push_back(entry.second); + + // normalize fAbundance (in mass only) if anyAbundanceGivenInMol + if (anyAbundanceGivenInMol) { + double sumMass = 0; + for (auto& nucl : fNuclei) sumMass += nucl.fAbundance; + for (auto& nucl : fNuclei) nucl.fAbundance /= sumMass; + } + + for (auto& nucl : fNuclei) nucl.PrintNucleus(); } /////////////////////////////////////////////// /// \brief Get recoil spectra for a given WIMP /// mass and cross section -/// Better performance version +/// Better performance version (velocity integral +/// is done only once for all recoil energies). /// std::map TRestWimpSensitivity::GetRecoilSpectra(const double wimpMass, const double crossSection) { diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx index 816c95b..4fe6cd6 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -207,3 +207,87 @@ const double TRestWimpUtils::GetQuenchingFactor(const double recoilEnergy, const return (1 + g) / g; } + +////////////////////////////////////////////////// +/// \brief Parse a chemical compound into a map +/// of elements and coefficients. The compound +/// can contain parentheses to indicate a +/// subCompound. The subCompound can have a +/// coefficient after the closing parenthesis. +/// The function is recursive, so subCompounds +/// can contain subCompounds. Examples: +/// "H2O" -> {"H": 2, "O": 1} +/// "Ne98(C4H10)2" -> {"Ne": 98, "C": 8, "H": 20} +/// "C6H5CH(CH3)2" -> {"C": 9, "H": 12} +/// +std::map TRestWimpUtils::ParseChemicalCompound(const std::string& compound) { + std::map elementMap; + std::string elementName; + // Get the number of opnening and closing parentheses + std::pair parenthesisCount = {0, 0}; + for (size_t i = 0; i < compound.size();) { + if (compound[i] == '(') parenthesisCount.first++; + if (compound[i] == ')') parenthesisCount.second++; + i++; + } + if (parenthesisCount.first != parenthesisCount.second) { + std::cout << "Error: Parentheses in compound " << compound << " do not match." << std::endl; + return elementMap; // empty map + } + + for (size_t i = 0; i < compound.size();) { + int coefficient = 1; + // Check for uppercase letter (start of an element) + if (std::isupper(compound[i])) { + elementName = compound[i]; + i++; + + // Check for lowercase letters (additional characters in element name) + while (i < compound.size() && std::islower(compound[i])) { + elementName += compound[i]; + i++; + } + + // Check for a number (coefficient) + if (i < compound.size() && std::isdigit(compound[i])) { + coefficient = 0; + while (i < compound.size() && std::isdigit(compound[i])) { + coefficient = coefficient * 10 + (compound[i] - '0'); + i++; + } + } + // Add the element and coefficient to the map + elementMap[elementName] += coefficient; + } else if (compound[i] == '(') { // Check for a subCompound inside parentheses + i++; + std::string subCompound; + while (i < compound.size()) { + if (compound[i] == ')') { + if (parenthesisCount.second > 1) + parenthesisCount.second--; + else + break; + } + subCompound += compound[i]; + i++; + } + i++; // Move past the closing parenthesis + // Find the subscript after the closing parenthesis + coefficient = 1; + if (i < compound.size() && std::isdigit(compound[i])) { + coefficient = 0; + while (i < compound.size() && std::isdigit(compound[i])) { + coefficient = coefficient * 10 + (compound[i] - '0'); + i++; + } + } + // Recursively call the function to handle the subCompound + std::map subElementMap = ParseChemicalCompound(subCompound); + for (auto& pair : subElementMap) { + elementMap[pair.first] += pair.second * coefficient; + } + } else + i++; + } + return elementMap; +}