Merge pull request #12758 from mantidproject/12741_absorption_coefficients

Fix normalization of atoms for compounds

Code/Mantid/Framework/DataHandling/src/SetSampleMaterial.cpp

@@ -57,7 +57,7 @@ void SetSampleMaterial::init() {
                   "formulas per cubic angstrom will be used instead of "
                   "calculated");
   declareProperty("ZParameter", EMPTY_DBL(), mustBePositive,
-                  "Number of atoms in the unit cell");
+                  "Number of formula units in unit cell");
   declareProperty("UnitCellVolume", EMPTY_DBL(), mustBePositive,
                   "Unit cell volume in Angstoms^3. Will be calculated from the "
                   "OrientedLattice if not supplied.");
@@ -202,20 +202,6 @@ void SetSampleMaterial::exec() {
   // determine the sample number density
   double rho = getProperty("SampleNumberDensity"); // in Angstroms-3
   double zParameter = getProperty("ZParameter");   // number of atoms
-  if (isEmpty(rho)) {
-    double unitCellVolume = getProperty("UnitCellVolume"); // in Angstroms^3
-
-    // get the unit cell volume from the workspace if it isn't set
-    if (isEmpty(unitCellVolume) && expInfo->sample().hasOrientedLattice()) {
-      unitCellVolume = expInfo->sample().getOrientedLattice().volume();
-      g_log.notice() << "found unit cell volume " << unitCellVolume
-                     << " Angstrom^-3\n";
-    }
-    // density is just number of atoms in the unit cell
-    // ...but only calculate it if you have both numbers
-    if ((!isEmpty(zParameter)) && (!isEmpty(unitCellVolume)))
-      rho = zParameter / unitCellVolume;
-  }
 
   // get the scattering information - this will override table values
   double coh_xs = getProperty("CoherentXSection");         // in barns
@@ -273,6 +259,21 @@ void SetSampleMaterial::exec() {
       // normalize the accumulated number by the number of atoms
       neutron = (1. / numAtoms) *
                 neutron; // funny syntax b/c of operators in neutron atom
+      if (isEmpty(rho)) {
+        double unitCellVolume = getProperty("UnitCellVolume"); // in Angstroms^3
+
+        // get the unit cell volume from the workspace if it isn't set
+        if (isEmpty(unitCellVolume) && expInfo->sample().hasOrientedLattice()) {
+          unitCellVolume = expInfo->sample().getOrientedLattice().volume();
+          g_log.notice() << "found unit cell volume " << unitCellVolume
+                         << " Angstrom^-3\n";
+        }
+        // density is just number of atoms in the unit cell
+        // ...but only calculate it if you have both numbers
+        if ((!isEmpty(zParameter)) && (!isEmpty(unitCellVolume)))
+          rho = numAtoms * zParameter / unitCellVolume;
+      }
+
       b_avg = b_avg / numAtoms;
       b_sq_avg = b_sq_avg / numAtoms;
 
@@ -289,6 +290,20 @@ void SetSampleMaterial::exec() {
     fixNeutron(neutron, coh_xs, inc_xs, sigma_atten, sigma_s);
 
     // create the material
+    if (isEmpty(rho)) {
+      double unitCellVolume = getProperty("UnitCellVolume"); // in Angstroms^3
+
+      // get the unit cell volume from the workspace if it isn't set
+      if (isEmpty(unitCellVolume) && expInfo->sample().hasOrientedLattice()) {
+        unitCellVolume = expInfo->sample().getOrientedLattice().volume();
+        g_log.notice() << "found unit cell volume " << unitCellVolume
+                       << " Angstrom^-3\n";
+      }
+      // density is just number of atoms in the unit cell
+      // ...but only calculate it if you have both numbers
+      if ((!isEmpty(zParameter)) && (!isEmpty(unitCellVolume)))
+        rho = zParameter / unitCellVolume;
+    }
     mat.reset(new Material(chemicalSymbol, neutron, rho));
   }
 

Code/Mantid/Framework/DataHandling/test/SetSampleMaterialTest.h

@@ -54,9 +54,9 @@ public:
 
 	  TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("InputWorkspace", wsName) );
     TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("ChemicalFormula","Al2-O3") );
-    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("SampleNumberDensity","0.0236649") );
-    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("ScatteringXSection","15.7048") );
-    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("AttenuationXSection","0.46257") );
+    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("SampleNumberDensity","0.1183245") );
+    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("ScatteringXSection","3.1404") );
+    TS_ASSERT_THROWS_NOTHING( setmat->setPropertyValue("AttenuationXSection","0.0925") );
     TS_ASSERT_THROWS_NOTHING( setmat->execute() );
     TS_ASSERT( setmat->isExecuted() );
 
@@ -70,9 +70,9 @@ public:
 
 	//can get away with holding pointer as it is an inout ws property
     const Material *m_sampleMaterial = &(testWS->sample().getMaterial());
-    TS_ASSERT_DELTA( m_sampleMaterial->numberDensity(), 0.0236649, 0.0001 );
-    TS_ASSERT_DELTA( m_sampleMaterial->totalScatterXSection(NeutronAtom::ReferenceLambda), 15.7048, 0.0001);
-    TS_ASSERT_DELTA( m_sampleMaterial->absorbXSection(NeutronAtom::ReferenceLambda), 0.46257, 0.0001);
+    TS_ASSERT_DELTA( m_sampleMaterial->numberDensity(), 0.1183245, 0.0001 );
+    TS_ASSERT_DELTA( m_sampleMaterial->totalScatterXSection(NeutronAtom::ReferenceLambda), 3.1404, 0.0001);
+    TS_ASSERT_DELTA( m_sampleMaterial->absorbXSection(NeutronAtom::ReferenceLambda), 0.0925, 0.0001);
 
     checkOutputProperties(setmat,m_sampleMaterial);
 
@@ -102,7 +102,7 @@ public:
     TS_ASSERT( setmat->isExecuted() );
 
     const Material *m_sampleMaterial = &(testWS->sample().getMaterial());
-    TS_ASSERT_DELTA( m_sampleMaterial->numberDensity(), 0.0236649, 0.0001 );
+    TS_ASSERT_DELTA( m_sampleMaterial->numberDensity(), 0.1183245, 0.0001 );
     TS_ASSERT_DELTA( m_sampleMaterial->totalScatterXSection(NeutronAtom::ReferenceLambda), 3.1404, 0.0001);
     TS_ASSERT_DELTA( m_sampleMaterial->absorbXSection(NeutronAtom::ReferenceLambda), 0.0925, 0.0001);
 	  

Code/Mantid/docs/source/algorithms/SetSampleMaterial-v1.rst

@@ -61,7 +61,7 @@ Number Density
 
 The number density is defined as
 
-.. math:: \rho_n = \frac{ZParameter}{UnitCellVolume}
+.. math:: \rho_n = \frac{N_{atoms}ZParameter}{UnitCellVolume}
 
 It can can be generated in one of two ways:
 
@@ -69,5 +69,22 @@ It can can be generated in one of two ways:
 2. Specifying the ``ZParameter`` and the ``UnitCellVolume`` (or letting
    the algorithm calculate it from the OrientedLattice on the 
    ``InputWorkspace``).
+   
+Linear Absorption Coefficients
+##############################
 
+.. math:: \mu_s = \rho_n \frac{1}{N_{atoms}}\sum_{i}s_{i}n_{i} \text{ units of 1/cm}
+.. math:: s = \sigma_{total scattering}
+.. math:: \mu_a = \rho_n \frac{1}{N_{atoms}}\sum_{i}a_{i}n_{i} \text{ units of 1/cm}
+.. math:: a = \sigma_{absorption} (\lambda=1.8)
+
+References
+----------
+
+The data used in this algorithm comes from the following paper.
+
+#. Varley F. Sears, *Neutron scattering lengths and cross sections*, Neutron News **3:3** (1992) 26
+   `doi: 10.1080/10448639208218770`_
+
+      
 .. categories::