Merge pull request #11 from rest-for-physics/alvaroezq_Volume

Functionality for setting nuclei targets abundances in mol (volume) and as chemical compound

examples/WIMP.rml

@@ -10,7 +10,7 @@
         <parameter name="escapeVelocity" value="544"/>
         <parameter name="exposure" value="116.8"/>
         <parameter name="background" value="1"/>
-        <parameter name="energySpectra" value="(0,2)"/>
+        <parameter name="energySpectra" value="(0,4)"/>
         <parameter name="energySpectraStep" value="0.01"/>
         <parameter name="energyRange" value="(0.1,1.1)"/>
         <parameter name="useQuenchingFactor" value="true"/>

examples/WIMP_compound_1.rml

@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<wimp>
+    <TRestWimpSensitivity name="WIMPSensitivity" title="WIMP Sensitivity" verboseLevel="info">
+        <addElement nucleusName="Ar" anum="39.948" znum="18" abundanceInMol="0.99"/>
+        <addCompound compoundName="C4H10" abundanceInMol="0.01">
+            <addElement nucleusName="C" anum="12.0107" znum="6" />
+            <addElement nucleusName="H" anum="1.00784" znum="1" />
+        </addCompound>
+        <parameter name="wimpDensity" value="0.3"/>
+        <parameter name="labVelocity" value="232"/>
+        <parameter name="rmsVelocity" value="220"/>
+        <parameter name="escapeVelocity" value="544"/>
+        <parameter name="exposure" value="116.8"/>
+        <parameter name="background" value="100"/>
+        <parameter name="energySpectra" value="(0,8)"/>
+        <parameter name="energySpectraStep" value="0.01"/>
+        <parameter name="energyRange" value="(0.05,1.05)"/>
+        <parameter name="useQuenchingFactor" value="true"/>
+    </TRestWimpSensitivity>
+</wimp>

examples/WIMP_compound_2.rml

@@ -0,0 +1,20 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<wimp>
+    <TRestWimpSensitivity name="WIMPSensitivity" title="WIMP Sensitivity" verboseLevel="info">
+        <addCompound compoundName="Ar99(C4H10)1" abundanceInMol="1">
+            <addElement nucleusName="Ar" anum="39.948" znum="18"/>
+            <addElement nucleusName="C" anum="12.0107" znum="6" />
+            <addElement nucleusName="H" anum="1.00784" znum="1" />
+        </addCompound>
+        <parameter name="wimpDensity" value="0.3"/>
+        <parameter name="labVelocity" value="232"/>
+        <parameter name="rmsVelocity" value="220"/>
+        <parameter name="escapeVelocity" value="544"/>
+        <parameter name="exposure" value="116.8"/>
+        <parameter name="background" value="100"/>
+        <parameter name="energySpectra" value="(0,8)"/>
+        <parameter name="energySpectraStep" value="0.01"/>
+        <parameter name="energyRange" value="(0.05,1.05)"/>
+        <parameter name="useQuenchingFactor" value="true"/>
+    </TRestWimpSensitivity>
+</wimp>

inc/TRestWimpNucleus.h

@@ -38,16 +38,19 @@ 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();
 
     // Destructor
     virtual ~TRestWimpNucleus();
 
-    ClassDef(TRestWimpNucleus, 1);
+    ClassDef(TRestWimpNucleus, 2);
 };
 
 #endif

inc/TRestWimpUtils.h

@@ -21,6 +21,7 @@
  *************************************************************************/
 
 #include <iostream>
+#include <map>
 
 #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<std::string, int> ParseChemicalCompound(const std::string& compound);
 }  // namespace TRestWimpUtils
 
 #endif

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;
+}

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 @@
 ///        <parameter name="useQuenchingFactor" value="true"/>
 ///    </TRestWimpSensitivity>
 /// \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<std::tuple<TString, Double_t, Int_t>, 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<std::string, TH1D*> TRestWimpSensitivity::GetRecoilSpectra(const double wimpMass,
                                                                     const double crossSection) {

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<std::string, int> TRestWimpUtils::ParseChemicalCompound(const std::string& compound) {
+    std::map<std::string, int> elementMap;
+    std::string elementName;
+    // Get the number of opnening and closing parentheses
+    std::pair<int, int> 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<std::string, int> subElementMap = ParseChemicalCompound(subCompound);
+            for (auto& pair : subElementMap) {
+                elementMap[pair.first] += pair.second * coefficient;
+            }
+        } else
+            i++;
+    }
+    return elementMap;
+}