Index: src/sas/sasgui/perspectives/fitting/fitting.py
===================================================================
--- src/sas/sasgui/perspectives/fitting/fitting.py	(revision 7432acb919f4013119ca918a73350297e0d4c93f)
+++ src/sas/sasgui/perspectives/fitting/fitting.py	(revision 24b382118c3145d493e1f76e2b308e6aa6f0759c)
@@ -359,5 +359,5 @@
                                    'Add a new model function')
         wx.EVT_MENU(owner, wx_id, self.make_new_model)
-        
+
         wx_id = wx.NewId()
         self.edit_model_menu.Append(wx_id, 'Sum|Multi(p1, p2)',
@@ -381,5 +381,5 @@
           '(Re)Load all models present in user plugin_models folder')
         wx.EVT_MENU(owner, wx_id, self.load_plugin_models)
-                
+   
     def set_edit_menu_helper(self, owner=None, menu=None):
         """
@@ -1748,11 +1748,11 @@
                                           dy=unsmeared_error)
             # Comment this out until we can get P*S models with correctly populated parameters
-            #if sq_model is not None and pq_model is not None:
-            #    self.create_theory_1D(x, sq_model, page_id, model, data, state,
-            #                          data_description=model.name + " S(q)",
-            #                          data_id=str(page_id) + " " + data.name + " S(q)")
-            #    self.create_theory_1D(x, pq_model, page_id, model, data, state,
-            #                          data_description=model.name + " P(q)",
-            #                          data_id=str(page_id) + " " + data.name + " P(q)")
+            if sq_model is not None and pq_model is not None:
+                self.create_theory_1D(x, sq_model, page_id, model, data, state,
+                                      data_description=model.name + " S(q)",
+                                      data_id=str(page_id) + " " + data.name + " S(q)")
+                self.create_theory_1D(x, pq_model, page_id, model, data, state,
+                                      data_description=model.name + " P(q)",
+                                      data_id=str(page_id) + " " + data.name + " P(q)")
 
             current_pg = self.fit_panel.get_page_by_id(page_id)
@@ -1910,5 +1910,5 @@
                 ## and may be the cause of other noted instabilities
                 ##
-                ##    -PDB August 12, 2014 
+                ##    -PDB August 12, 2014
                 while self.calc_2D.isrunning():
                     time.sleep(0.1)
@@ -1952,5 +1952,5 @@
             if (self.calc_1D is not None) and self.calc_1D.isrunning():
                 self.calc_1D.stop()
-                ## stop just raises the flag -- the thread is supposed to 
+                ## stop just raises the flag -- the thread is supposed to
                 ## then kill itself but cannot.  Paul Kienzle came up with
                 ## this fix to prevent threads from stepping on each other
@@ -1964,5 +1964,5 @@
                 ## a request to stop the computation.
                 ## It seems thus that the whole thread approach used here
-                ## May need rethinking  
+                ## May need rethinking
                 ##
                 ##    -PDB August 12, 2014
@@ -2129,5 +2129,5 @@
             residuals.dxw = None
             residuals.ytransform = 'y'
-            # For latter scale changes 
+            # For latter scale changes
             residuals.xaxis('\\rm{Q} ', 'A^{-1}')
             residuals.yaxis('\\rm{Residuals} ', 'normalized')
Index: src/sas/sasgui/perspectives/fitting/model_thread.py
===================================================================
--- src/sas/sasgui/perspectives/fitting/model_thread.py	(revision 7432acb919f4013119ca918a73350297e0d4c93f)
+++ src/sas/sasgui/perspectives/fitting/model_thread.py	(revision 24b382118c3145d493e1f76e2b308e6aa6f0759c)
@@ -71,5 +71,5 @@
                     (self.data.qy_data * self.data.qy_data))
 
-        # For theory, qmax is based on 1d qmax 
+        # For theory, qmax is based on 1d qmax
         # so that must be mulitified by sqrt(2) to get actual max for 2d
         index_model = (self.qmin <= radius) & (radius <= self.qmax)
@@ -91,9 +91,10 @@
                 self.data.qy_data[index_model]
             ])
-        output = np.zeros(len(self.data.qx_data))
+        # Initialize output to NaN so masked elements do not get plotted.
+        output = np.empty_like(self.data.qx_data)
         # output default is None
         # This method is to distinguish between masked
         #point(nan) and data point = 0.
-        output = output / output
+        output[:] = np.NaN
         # set value for self.mask==True, else still None to Plottools
         output[index_model] = value
@@ -198,23 +199,27 @@
             output[index] = self.model.evalDistribution(self.data.x[index])
 
+        x=self.data.x[index]
+        y=output[index]
         sq_values = None
         pq_values = None
-        s_model = None
-        p_model = None
         if isinstance(self.model, MultiplicationModel):
             s_model = self.model.s_model
             p_model = self.model.p_model
-        elif hasattr(self.model, "get_composition_models"):
-            p_model, s_model = self.model.get_composition_models()
-
-        if p_model is not None and s_model is not None:
-            sq_values = np.zeros((len(self.data.x)))
-            pq_values = np.zeros((len(self.data.x)))
-            sq_values[index] = s_model.evalDistribution(self.data.x[index])
-            pq_values[index] = p_model.evalDistribution(self.data.x[index])
+            #sq_values = np.zeros((len(self.data.x)))
+            #pq_values = np.zeros((len(self.data.x)))
+            #sq_values[index] = s_model.evalDistribution(self.data.x[index])
+            #pq_values[index] = p_model.evalDistribution(self.data.x[index])
+            sq_values = s_model.evalDistribution(x)
+            pq_values = p_model.evalDistribution(x)
+        elif hasattr(self.model, "calc_composition_models"):
+            results = self.model.calc_composition_models(x)
+            if results is not None:
+                sq_values = results[0]
+                pq_values = results[1]
+
 
         elapsed = time.time() - self.starttime
 
-        self.complete(x=self.data.x[index], y=output[index],
+        self.complete(x=x, y=y,
                       page_id=self.page_id,
                       state=self.state,
