diff --git a/src/main.py b/src/main.py
index 369c53b..0ccfc30 100644
--- a/src/main.py
+++ b/src/main.py
@@ -8,7 +8,15 @@ if __name__ == "__main__":
     parser.add_argument('--num_runs', type=int, default=10, help='number of simulations to run with a fixed set of parameters')
     parser.add_argument('--min_runs', type=int, default=5, help='minimum number of successful runs needed to calculate statistics')
     parser.add_argument('--confidence_level', type=float, default=0.95, help='confidence level')
+    parser.add_argument('--best_c', action=argparse.BooleanOptionalAction, help='only run simulations for lambda=1 and determine the best C')
 
     args = parser.parse_args()
-    simulation_wrapper(args.simulation_time, args.num_runs, args.min_runs, args.confidence_level)
+
+    if args.best_c:
+        lambda_vals = [1,]
+    else:
+        lambda_vals = [l/100 for l in range(1, 301)]  # λ from 0.01 to 3.00
+
+
+    simulation_wrapper(args.simulation_time, args.num_runs, args.min_runs, args.confidence_level, lambda_vals)
 
diff --git a/src/simulation.py b/src/simulation.py
index ed2f221..dc5eb3b 100644
--- a/src/simulation.py
+++ b/src/simulation.py
@@ -149,13 +149,13 @@ def run_single_simulation(args):
     except ValueError:  # Loss rate too high
         return None
 
-def simulation_wrapper(simulation_time, num_runs, min_runs, confidence_level):
+def simulation_wrapper(simulation_time, num_runs, min_runs, confidence_level, lambda_vals):
     C_values = [1, 2, 3, 6]
 
-    lambda_vals = [l/100 for l in range(1, 301)]  # λ from 0.01 to 3.00
-
     plt.figure(figsize=(12, 8))
 
+    mean_at_lambda1 = {}
+
     with Pool() as pool: # pool of workers
         for c in C_values:
             lambda_points = []
@@ -202,13 +202,31 @@ def simulation_wrapper(simulation_time, num_runs, min_runs, confidence_level):
                     print(f"Stopped at λ={lambda_val:.2f} - no successful run")
                     break
 
-
             plt.plot(lambda_points, means, label=f'C={c}')
             plt.fill_between(lambda_points, ci_lower, ci_upper, alpha=0.2)
 
-    plt.xlabel('Arrival Rate (λ)')
-    plt.ylabel('Mean Response Time')
-    plt.title(f'Mean Response Time vs Arrival Rate ({num_runs} runs, 95% CI)')
-    plt.legend()
-    plt.grid(True)
-    plt.show()
+            # store response time for lamba = 1
+            if 1 in lambda_points:
+                idx = lambda_points.index(1)
+                mean_at_lambda1[c] = (means[idx], ci_lower[idx], ci_upper[idx])
+
+    # determine optimal C for lamba = 1
+    if mean_at_lambda1:
+        sorted_C = sorted(mean_at_lambda1.items(), key=lambda item: item[1][0])
+        (best_C, (mean1, lower1, upper1)), (_, (_, lower2, _)) = sorted_C[:2]
+        if upper1 < lower2:
+            print(f"Optimal C at λ=1 is {best_C} with Mean RT = {mean1:.2f} (non-overlapping CIs)")
+        else:
+            print("Confidence intervals overlap between the two best C.")
+    else:
+        print("\nNo valid λ=1 data for any C.")
+
+
+    # plot curves
+    if len(lambda_vals)>1:
+        plt.xlabel('Arrival Rate (λ)')
+        plt.ylabel('Mean Response Time')
+        plt.title(f'Mean Response Time vs Arrival Rate ({num_runs} runs, 95% CI)')
+        plt.legend()
+        plt.grid(True)
+        plt.show()