diff --git a/plotting.py b/plotting.py
index f5793624aaa240d32e1313445485d7e8b2fa9bd6..05110526959b5bd15e18bf1627bd68ba8026df4a 100644
--- a/plotting.py
+++ b/plotting.py
@@ -226,7 +226,7 @@ def plot_hist_2D(plotname, xedges, yedges, hist, varx_label=None, vary_label=Non
         extraopts.update(norm=matplotlib.colors.LogNorm(vmin=np.min(hist[hist>0]), vmax=np.max(hist)))
 
     ax.set_facecolor("black")
-    pcm = ax.pcolormesh(X, Y, hist, cmap="inferno", **extraopts)
+    pcm = ax.pcolormesh(X, Y, hist, cmap="inferno", linewidth=0, rasterized=True, **extraopts)
     cbar = fig.colorbar(pcm, ax=ax)
     cbar.set_label(zlabel)
     ax.set_ylabel(vary_label)
@@ -254,12 +254,13 @@ def plot_cond_avg_actmax_2D(plotname, model, layer, neuron, ranges,
                             varx_index, vary_index,
                             nbinsx, xmin, xmax, nbinsy, ymin, ymax,
                             scaler=None,
+                            ntries=20,
                             **kwargs):
 
     xedges = np.linspace(xmin, xmax, nbinsx)
     yedges = np.linspace(ymin, ymax, nbinsy)
 
-    hist = np.zeros(nbinsx*nbinsy).reshape(nbinsx, nbinsy)
+    hist = np.zeros(int(nbinsx*nbinsy)).reshape(int(nbinsx), int(nbinsy))
 
     gradient_function = get_grad_function(model, layer, neuron)
 
@@ -272,9 +273,11 @@ def plot_cond_avg_actmax_2D(plotname, model, layer, neuron, ranges,
                 random_event[0][index] = val
             if scaler is not None:
                 random_event = scaler.transform(random_event)
-            act = np.mean([max_activation_wrt_input(gradient_function, random_event, maxit=1, const_indices=[varx_index, vary_index])[0][0] for i in range(20)])
+            act = np.mean([max_activation_wrt_input(gradient_function, random_event, maxit=1, const_indices=[varx_index, vary_index])[0][0] for i in range(ntries)])
             hist[ix][iy] = act
 
+    hist = hist.T
+
     plot_hist_2D(plotname, xedges, yedges, hist, zlabel="Neuron output", **kwargs)
 
 
diff --git a/scripts/plot_single_neuron.py b/scripts/plot_single_neuron.py
index 6f128199d8b11c4d4b86b1ce4aabf3127593b456..ce6728f62a3159ed12579aeef46efd31c5131718 100755
--- a/scripts/plot_single_neuron.py
+++ b/scripts/plot_single_neuron.py
@@ -10,7 +10,8 @@ from KerasROOTClassification.plotting import (
     get_mean_event,
     plot_NN_vs_var_2D,
     plot_profile_2D,
-    plot_hist_2D_events
+    plot_hist_2D_events,
+    plot_cond_avg_actmax_2D
 )
 from KerasROOTClassification.tfhelpers import get_single_neuron_function
 
@@ -20,7 +21,7 @@ parser.add_argument("output_filename")
 parser.add_argument("varx")
 parser.add_argument("vary")
 parser.add_argument("-m", "--mode",
-                    choices=["mean_sig", "mean_bkg", "profile_sig", "profile_bkg", "hist_sig", "hist_bkg"],
+                    choices=["mean_sig", "mean_bkg", "profile_sig", "profile_bkg", "hist_sig", "hist_bkg", "cond_actmax"],
                     default="mean_sig")
 parser.add_argument("-l", "--layer", type=int, help="Layer index (takes last layer by default)")
 parser.add_argument("-n", "--neuron", type=int, default=0, help="Neuron index (takes first neuron by default)")
@@ -29,7 +30,8 @@ parser.add_argument("--contour", action="store_true", help="Interpolate with con
 parser.add_argument("-b", "--nbins", default=20, type=int, help="Number of bins in x and y direction")
 parser.add_argument("-x", "--xrange", type=float, nargs="+", help="xrange (low, high)")
 parser.add_argument("-y", "--yrange", type=float, nargs="+", help="yrange (low, high)")
-parser.add_argument("-p", "--profile-metric", help="metric for profile modes", default="mean", choices=["mean", "average", "max"])
+parser.add_argument("-p", "--profile-metric", help="metric for profile modes", default="average", choices=["mean", "average", "max"])
+parser.add_argument("--ntries-cond-actmax", help="number of random events to be maximised and averaged per bin", default=20, type=int)
 
 args = parser.parse_args()
 
@@ -111,6 +113,7 @@ elif args.mode.startswith("profile"):
         ymin=vary_range[0], ymax=vary_range[1], nbinsy=vary_range[2],
         metric=metric_dict[args.profile_metric],
         varx_label=varx_label, vary_label=vary_label,
+        log=args.log,
         **opt_kwargs
     )
 
@@ -134,3 +137,23 @@ elif args.mode.startswith("hist"):
         varx_label=varx_label, vary_label=vary_label,
         log=args.log,
     )
+
+elif args.mode.startswith("cond_actmax"):
+
+    x_test_scaled = c.scaler.transform(c.x_test)
+
+    # ranges in which to sample the random events
+    ranges = [np.percentile(x_test_scaled[:,var_index], [1,99]) for var_index in range(len(c.branches))]
+
+    plot_cond_avg_actmax_2D(
+        args.output_filename,
+        c.model, layer, neuron, ranges,
+        varx_index,
+        vary_index,
+        xmin=varx_range[0], xmax=varx_range[1], nbinsx=varx_range[2],
+        ymin=vary_range[0], ymax=vary_range[1], nbinsy=vary_range[2],
+        scaler=c.scaler,
+        ntries=args.ntries_cond_actmax,
+        varx_label=varx_label, vary_label=vary_label,
+        log=args.log,
+    )