From e53037c8b0057234dfe5ff6ed84dacbaa417c7f4 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Thu, 2 Nov 2023 12:58:37 +0100 Subject: [PATCH 1/6] Adding functionality for defining target material with abundance in mol/volume and chemical compounds --- examples/WIMP_compound_1.rml | 20 +++++++ examples/WIMP_compound_2.rml | 20 +++++++ inc/TRestWimpNucleus.h | 5 +- inc/TRestWimpUtils.h | 3 +- src/TRestWimpNucleus.cxx | 27 ++++++++- src/TRestWimpSensitivity.cxx | 104 ++++++++++++++++++++++++++++++++++- src/TRestWimpUtils.cxx | 66 ++++++++++++++++++++++ 7 files changed, 241 insertions(+), 4 deletions(-) create mode 100644 examples/WIMP_compound_1.rml create mode 100644 examples/WIMP_compound_2.rml 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..0e0bec6 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) percentage + 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..d578618 100644 --- a/src/TRestWimpNucleus.cxx +++ b/src/TRestWimpNucleus.cxx @@ -40,12 +40,19 @@ /// #include "TRestWimpNucleus.h" +#include "TRestWimpUtils.h" #include "TRestMetadata.h" ClassImp(TRestWimpNucleus); -TRestWimpNucleus::TRestWimpNucleus() {} +TRestWimpNucleus::TRestWimpNucleus() { + fNucleusName = ""; + fAnum = 0; + fZnum = 0; + fAbundance = 0; + fAbundanceMol = 0; +} TRestWimpNucleus::~TRestWimpNucleus() {} @@ -55,5 +62,23 @@ 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; } + +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; +} \ No newline at end of file diff --git a/src/TRestWimpSensitivity.cxx b/src/TRestWimpSensitivity.cxx index 7987e6c..4e0aa4a 100644 --- a/src/TRestWimpSensitivity.cxx +++ b/src/TRestWimpSensitivity.cxx @@ -68,6 +68,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,17 +167,113 @@ 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(); } /////////////////////////////////////////////// diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx index 816c95b..43037b0 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -207,3 +207,69 @@ 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. +/// +std::map TRestWimpUtils::ParseChemicalCompound(const std::string& compound) { + std::map elementMap; + std::string elementName; + int coefficient = 1; + + for (size_t i = 0; i < compound.size(); ) { + // 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() && compound[i] != ')') { + 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; +} \ No newline at end of file From 0bdf8a792c94459980eee255cdab304ecf2d9e79 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Mon, 6 Nov 2023 11:47:45 +0100 Subject: [PATCH 2/6] pre-commit fixes --- src/TRestWimpNucleus.cxx | 7 +++---- src/TRestWimpSensitivity.cxx | 33 ++++++++++++++++++--------------- src/TRestWimpUtils.cxx | 14 ++++++-------- 3 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/TRestWimpNucleus.cxx b/src/TRestWimpNucleus.cxx index d578618..002027a 100644 --- a/src/TRestWimpNucleus.cxx +++ b/src/TRestWimpNucleus.cxx @@ -40,9 +40,9 @@ /// #include "TRestWimpNucleus.h" -#include "TRestWimpUtils.h" #include "TRestMetadata.h" +#include "TRestWimpUtils.h" ClassImp(TRestWimpNucleus); @@ -67,7 +67,6 @@ void TRestWimpNucleus::PrintNucleus() { } int TRestWimpNucleus::GetStechiometricFactorFromCompound(const std::string& compound) { - auto elementMap = TRestWimpUtils::ParseChemicalCompound(compound); int stechiometricFactor = 0; @@ -79,6 +78,6 @@ int TRestWimpNucleus::GetStechiometricFactorFromCompound(const std::string& comp } if (stechiometricFactor == 0) RESTWarning << "No nucleus " << fNucleusName.Data() << " founnd in compound " << compound << RESTendl; - + return stechiometricFactor; -} \ No newline at end of file +} diff --git a/src/TRestWimpSensitivity.cxx b/src/TRestWimpSensitivity.cxx index 4e0aa4a..467f346 100644 --- a/src/TRestWimpSensitivity.cxx +++ b/src/TRestWimpSensitivity.cxx @@ -177,10 +177,12 @@ void TRestWimpSensitivity::ReadNuclei() { nucleus.fAnum = StringToDouble(GetFieldValue("anum", ElementDef)); nucleus.fZnum = StringToInteger(GetFieldValue("znum", ElementDef)); - std::string el = !ElementDef->Attribute("abundance") ? "Not defined" : ElementDef->Attribute("abundance"); + 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"); + el = !ElementDef->Attribute("abundanceInMol") ? "Not defined" + : ElementDef->Attribute("abundanceInMol"); if (!(el.empty() || el == "Not defined")) { nucleus.fAbundanceMol = StringToDouble(el); anyAbundanceGivenInMol = true; @@ -192,8 +194,8 @@ void TRestWimpSensitivity::ReadNuclei() { else if (nucleus.fAbundanceMol == 0) nucleus.fAbundanceMol = nucleus.fAbundance / nucleus.fAnum; else - RESTError << "abundance or abundanceInMol not defined for nucleus " - << nucleus.fNucleusName << RESTendl; + RESTError << "abundance or abundanceInMol not defined for nucleus " << nucleus.fNucleusName + << RESTendl; } fNuclei.emplace_back(nucleus); ElementDef = GetNextElement(ElementDef); @@ -202,16 +204,17 @@ void TRestWimpSensitivity::ReadNuclei() { // 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"); + 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"); + el = !CompoundDef->Attribute("abundanceInMol") ? "Not defined" + : CompoundDef->Attribute("abundanceInMol"); if (!(el.empty() || el == "Not defined")) { compoundAbundanceInMol = StringToDouble(el); compoundAbundanceGivenInMol = true; @@ -224,8 +227,8 @@ void TRestWimpSensitivity::ReadNuclei() { compoundAbundanceInMol = 1; } - // Read nuclei (inside compound) elements - TiXmlElement *ElementDef = GetElement("addElement", CompoundDef); + // Read nuclei (inside compound) elements + TiXmlElement* ElementDef = GetElement("addElement", CompoundDef); int i = 0; while (ElementDef) { i++; @@ -243,7 +246,7 @@ void TRestWimpSensitivity::ReadNuclei() { 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++) { + for (auto it = fNuclei.end() - i; it != fNuclei.end(); it++) { auto& nucleus = *it; int stechiometricFactor = nucleus.GetStechiometricFactorFromCompound(compoundName); nucleus.fAbundanceMol = compoundAbundanceInMol * stechiometricFactor; @@ -260,14 +263,14 @@ void TRestWimpSensitivity::ReadNuclei() { if (uniqueNucleiMap.find(key) != uniqueNucleiMap.end()) { uniqueNucleiMap[key].fAbundance += nucleus.fAbundance; uniqueNucleiMap[key].fAbundanceMol += nucleus.fAbundanceMol; - } else uniqueNucleiMap[key] = nucleus; + } else + uniqueNucleiMap[key] = nucleus; } fNuclei.clear(); - for (const auto& entry : uniqueNucleiMap) - fNuclei.push_back(entry.second); + for (const auto& entry : uniqueNucleiMap) fNuclei.push_back(entry.second); - //normalize fAbundance (in mass only) if anyAbundanceGivenInMol - if (anyAbundanceGivenInMol){ + // 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; diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx index 43037b0..a003181 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -221,8 +221,8 @@ std::map TRestWimpUtils::ParseChemicalCompound(const std::stri std::map elementMap; std::string elementName; int coefficient = 1; - - for (size_t i = 0; i < compound.size(); ) { + + for (size_t i = 0; i < compound.size();) { // Check for uppercase letter (start of an element) if (std::isupper(compound[i])) { elementName = compound[i]; @@ -244,15 +244,14 @@ std::map TRestWimpUtils::ParseChemicalCompound(const std::stri } // Add the element and coefficient to the map elementMap[elementName] += coefficient; - } - else if (compound[i] == '(') { // Check for a subCompound inside parentheses + } else if (compound[i] == '(') { // Check for a subCompound inside parentheses i++; std::string subCompound; while (i < compound.size() && compound[i] != ')') { subCompound += compound[i]; i++; } - i++; // Move past the closing parenthesis + i++; // Move past the closing parenthesis // Find the subscript after the closing parenthesis coefficient = 1; if (i < compound.size() && std::isdigit(compound[i])) { @@ -267,9 +266,8 @@ std::map TRestWimpUtils::ParseChemicalCompound(const std::stri for (auto& pair : subElementMap) { elementMap[pair.first] += pair.second * coefficient; } - } - else + } else i++; } return elementMap; -} \ No newline at end of file +} From d4e130f9b81a46fe02d771b76df7a6c84cc87863 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 8 Nov 2023 10:57:57 +0100 Subject: [PATCH 3/6] Improving documentation --- inc/TRestWimpNucleus.h | 2 +- src/TRestWimpNucleus.cxx | 6 ++++++ src/TRestWimpSensitivity.cxx | 5 ++++- src/TRestWimpUtils.cxx | 5 ++++- 4 files changed, 15 insertions(+), 3 deletions(-) diff --git a/inc/TRestWimpNucleus.h b/inc/TRestWimpNucleus.h index 0e0bec6..23f8a83 100644 --- a/inc/TRestWimpNucleus.h +++ b/inc/TRestWimpNucleus.h @@ -38,7 +38,7 @@ class TRestWimpNucleus { Int_t fZnum; /// Abundance, in mass percentage Double_t fAbundance; - /// Abundance, in mole (or volume) percentage + /// Abundance, in mole (or volume) Double_t fAbundanceMol; void PrintNucleus(); diff --git a/src/TRestWimpNucleus.cxx b/src/TRestWimpNucleus.cxx index 002027a..567c519 100644 --- a/src/TRestWimpNucleus.cxx +++ b/src/TRestWimpNucleus.cxx @@ -66,6 +66,12 @@ void TRestWimpNucleus::PrintNucleus() { 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); diff --git a/src/TRestWimpSensitivity.cxx b/src/TRestWimpSensitivity.cxx index 467f346..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 @@ -282,7 +284,8 @@ void TRestWimpSensitivity::ReadNuclei() { /////////////////////////////////////////////// /// \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 a003181..c7b5426 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -215,7 +215,10 @@ const double TRestWimpUtils::GetQuenchingFactor(const double recoilEnergy, const /// subCompound. The subCompound can have a /// coefficient after the closing parenthesis. /// The function is recursive, so subCompounds -/// can contain 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; From d3d3dd44112ec63f404505f9059336c217992e07 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 8 Nov 2023 11:04:10 +0100 Subject: [PATCH 4/6] Fixing bug TRestWimpUtils::ParseChemicalCompound --- src/TRestWimpUtils.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx index c7b5426..1f27b38 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -223,9 +223,9 @@ const double TRestWimpUtils::GetQuenchingFactor(const double recoilEnergy, const std::map TRestWimpUtils::ParseChemicalCompound(const std::string& compound) { std::map elementMap; std::string elementName; - int coefficient = 1; 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]; From de30ed988e36974ca17426667433c027dc0a0877 Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Wed, 8 Nov 2023 11:54:59 +0100 Subject: [PATCH 5/6] Fixing nested subcompound in ParseChemicalCompound() --- src/TRestWimpUtils.cxx | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/TRestWimpUtils.cxx b/src/TRestWimpUtils.cxx index 1f27b38..4fe6cd6 100644 --- a/src/TRestWimpUtils.cxx +++ b/src/TRestWimpUtils.cxx @@ -223,6 +223,17 @@ const double TRestWimpUtils::GetQuenchingFactor(const double recoilEnergy, const 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; @@ -250,7 +261,13 @@ std::map TRestWimpUtils::ParseChemicalCompound(const std::stri } else if (compound[i] == '(') { // Check for a subCompound inside parentheses i++; std::string subCompound; - while (i < compound.size() && compound[i] != ')') { + while (i < compound.size()) { + if (compound[i] == ')') { + if (parenthesisCount.second > 1) + parenthesisCount.second--; + else + break; + } subCompound += compound[i]; i++; } From 93039e32f1d6a486b236498e4320842b16c7cd9b Mon Sep 17 00:00:00 2001 From: Alvaro Ezquerro Date: Thu, 30 Nov 2023 11:36:18 +0100 Subject: [PATCH 6/6] Making energy spectra range wide enough --- examples/WIMP.rml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 @@ - +