diff --git a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08_qubit_spectroscopy.py b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08a_qubit_spectroscopy.py
similarity index 100%
rename from Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08_qubit_spectroscopy.py
rename to Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08a_qubit_spectroscopy.py
diff --git a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08b_qubit_spectroscopy_with_chirp.py b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08b_qubit_spectroscopy_with_chirp.py
new file mode 100644
index 000000000..0ab12529d
--- /dev/null
+++ b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/08b_qubit_spectroscopy_with_chirp.py
@@ -0,0 +1,245 @@
+"""
+        QUBIT SPECTROSCOPY w/ Chirp
+The goal of the script is to find the qubit transition by sweeping both the qubit pulse frequency and the magnetic field.
+The QUA program is divided into three sections:
+    1) step between the initialization point, idle point, and the measurement point using sticky elements (long timescale).
+    2) send the chirp pulse to drive the EDSR transition (short timescale).
+    3) measure the state of the qubit using either RF reflectometry or dc current sensing via PSB or Elzerman readout.
+A compensation pulse can be added to the long timescale sequence in order to ensure 0 DC voltage on the fast line of
+the bias-tee. Alternatively one can obtain the same result by changing the offset of the slow line of the bias-tee.
+
+In the current implementation, the magnetic field and LO frequency are being swept using the API of the relevant
+instruments in Python. For this reason the OPX program is paused at each iteration, the external parameters (B and f_LO)
+are updated in Python and then the QUA program is resumed to sweep the qubit intermediate frequency and measure the
+state of the dot.
+Also note that the qubit pulse is played at the end of the "idle" level whose duration is fixed.
+
+Note that providing a single magnetic field value will result in acquiring the 1D qubit spectroscopy at the specified
+B-field.
+
+Prerequisites:
+    - Readout calibration (resonance frequency for RF reflectometry and sensor operating point for DC current sensing).
+    - Setting the DC offsets of the external DC voltage source.
+    - Connecting the OPX to the fast line of the plunger gates.
+    - Having calibrated the initialization and readout point from the charge stability map and updated the configuration.
+
+Before proceeding to the next node:
+    - Identify the qubit frequency and update the configuration.
+"""
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool, wait_until_job_is_paused
+from qualang_tools.plot import interrupt_on_close
+from qualang_tools.loops import from_array
+import matplotlib.pyplot as plt
+from macros import RF_reflectometry_macro, DC_current_sensing_macro
+
+
+###################
+# The QUA program #
+###################
+
+n_avg = 100
+
+# Chirp parameters - Defined in configuration
+
+# The intermediate frequency sweep parameters
+f_min = 10 * u.MHz
+f_max = 251 * u.MHz
+df = 2000 * u.kHz
+IFs = np.arange(f_min, f_max + 0.1, df)
+# The LO frequency sweep parameters
+f_min_external = 4.501e9 - f_min
+f_max_external = 6.5e9 - f_max
+df_external = f_max - f_min
+lo_frequencies = np.arange(f_min_external, f_max_external + 0.1, df_external)
+# lo_frequencies = [6e9]
+# Total frequency vector
+frequencies = np.array(np.concatenate([IFs + lo_frequencies[i] for i in range(len(lo_frequencies))]))
+
+# Magnetic field in T
+# B_fields = np.arange(-5, 5, 0.1)
+B_fields = [0, 1, 2]
+
+# Delay in ns before stepping to the readout point after playing the qubit pulse - must be a multiple of 4ns and >= 16ns
+delay_before_readout = 16
+
+seq = OPX_virtual_gate_sequence(config, ["P1_sticky", "P2_sticky"])
+seq.add_points("initialization", level_init, duration_init)
+seq.add_points("idle", level_manip, duration_manip)
+seq.add_points("readout", level_readout, duration_readout)
+
+with program() as qubit_spectroscopy_prog:
+    n = declare(int)  # QUA integer used as an index for the averaging loop
+    f = declare(int)  # QUA variable for the qubit pulse duration
+    i = declare(int)  # QUA variable for the magnetic field sweep
+    j = declare(int)  # QUA variable for the lo frequency sweep
+    chirp_var = declare(int, value=chirp_rate)
+    n_st = declare_stream()  # Stream for the iteration number (progress bar)
+    with for_(i, 0, i < len(B_fields) + 1, i + 1):
+        with for_(j, 0, j < len(lo_frequencies), j + 1):
+            # pause() # Needs to be uncommented when not simulating
+            with for_(n, 0, n < n_avg, n + 1):  # The averaging loop
+                with for_(*from_array(f, IFs)):  # Loop over the qubit pulse amplitude
+                    update_frequency("qubit", f)
+
+                    # Navigate through the charge stability map
+                    seq.add_step(voltage_point_name="initialization")
+                    seq.add_step(
+                        voltage_point_name="idle", duration=chirp_duration + processing_time + delay_before_readout
+                    )  # Processing time is time it takes to calculate the chirp pulse
+                    seq.add_step(voltage_point_name="readout", duration=duration_readout)
+                    seq.add_compensation_pulse(duration=duration_compensation_pulse)
+
+                    # Drive the qubit by playing the MW pulse at the end of the manipulation step
+                    wait((duration_init) * u.ns, "qubit")  #
+                    play("chirp", "qubit", chirp=(chirp_var, chirp_units))
+
+                    # Measure the dot right after the qubit manipulation
+                    wait(
+                        (duration_init + chirp_duration + processing_time + delay_before_readout) * u.ns,
+                        "tank_circuit",
+                        "TIA",
+                    )  #
+                    I, Q, I_st, Q_st = RF_reflectometry_macro()
+                    dc_signal, dc_signal_st = DC_current_sensing_macro()
+
+                    seq.ramp_to_zero()
+        save(i, n_st)
+    # Stream processing section used to process the data before saving it.
+    with stream_processing():
+        n_st.save("iteration")
+        # Cast the data into a 2D matrix and performs a global averaging of the received 2D matrices together.
+        # RF reflectometry
+        I_st.buffer(len(IFs)).buffer(n_avg).map(FUNCTIONS.average()).buffer(len(lo_frequencies)).save_all("I")
+        Q_st.buffer(len(IFs)).buffer(n_avg).map(FUNCTIONS.average()).buffer(len(lo_frequencies)).save_all("Q")
+        # DC current sensing
+        dc_signal_st.buffer(len(IFs)).buffer(n_avg).map(FUNCTIONS.average()).buffer(len(lo_frequencies)).save_all(
+            "dc_signal"
+        )
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+simulate = True
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=10_000)  # In clock cycles = 4ns
+    # Simulate blocks python until the simulation is done
+    job = qmm.simulate(config, qubit_spectroscopy_prog, simulation_config)
+
+    # Get the waveform report
+    samples = job.get_simulated_samples()
+    waveform_report = job.get_simulated_waveform_report()
+    waveform_report.create_plot(samples, plot=True, save_path=None)
+
+    # Plot the simulated samples
+    plt.figure()
+    plt.subplot(211)
+    job.get_simulated_samples().con1.plot()
+    plt.axhline(level_init[0], color="k", linestyle="--")
+    plt.axhline(level_manip[0], color="k", linestyle="--")
+    plt.axhline(level_readout[0], color="k", linestyle="--")
+    plt.axhline(level_init[1], color="k", linestyle="--")
+    plt.axhline(level_manip[1], color="k", linestyle="--")
+    plt.axhline(level_readout[1], color="k", linestyle="--")
+    plt.yticks(
+        [
+            level_readout[1],
+            level_manip[1],
+            level_init[1],
+            0.0,
+            level_init[0],
+            level_manip[0],
+            level_readout[0],
+        ],
+        ["readout", "manip", "init", "0", "init", "manip", "readout"],
+    )
+    plt.legend("")
+    from macros import get_filtered_voltage
+
+    plt.subplot(212)
+    get_filtered_voltage(
+        job.get_simulated_samples().con1.analog["1"],
+        1e-9,
+        bias_tee_cut_off_frequency,
+        True,
+    )
+    plt.show()
+
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(qubit_spectroscopy_prog)
+    # Live plotting
+    fig = plt.figure()
+    interrupt_on_close(fig, job)  # Interrupts the job when closing the figure
+    for i in range(len(B_fields)):  # Loop over y-voltages
+        # TODO Update the magnetic field
+        for j in range(len(lo_frequencies)):
+            # TODO update the lo frequency
+            # Resume the QUA program (escape the 'pause' statement)
+            job.resume()
+            # Wait until the program reaches the 'pause' statement again, indicating that the QUA program is done
+            wait_until_job_is_paused(job)
+        if i == 0:
+            # Get results from QUA program and initialize live plotting
+            results = fetching_tool(job, data_list=["I", "Q", "dc_signal", "iteration"], mode="live")
+        # Fetch the data from the last OPX run corresponding to the current slow axis iteration
+        I, Q, DC_signal, iteration = results.fetch_all()
+        # Convert results into Volts
+        S = u.demod2volts(I + 1j * Q, reflectometry_readout_length)
+        R = np.abs(S)  # Amplitude
+        phase = np.angle(S)  # Phase
+        DC_signal = u.demod2volts(DC_signal, readout_len)
+        # Progress bar
+        progress_counter(iteration, len(B_fields))
+        # Plot data
+        if len(B_fields) > 1:
+            plt.subplot(121)
+            plt.cla()
+            plt.title(r"$R=\sqrt{I^2 + Q^2}$ [V]")
+            plt.pcolor(
+                frequencies / u.MHz,
+                B_fields[: iteration + 1],
+                np.reshape(R, (iteration + 1, len(frequencies))),
+            )
+            plt.xlabel("Qubit pulse frequency [MHz]")
+            plt.ylabel("B [mT]")
+            plt.subplot(122)
+            plt.cla()
+            plt.title("Phase [rad]")
+            plt.pcolor(
+                frequencies / u.MHz,
+                B_fields[: iteration + 1],
+                np.reshape(phase, (iteration + 1, len(frequencies))),
+            )
+            plt.xlabel("Qubit pulse frequency [MHz]")
+            plt.ylabel("B [mT]")
+            plt.tight_layout()
+            plt.pause(0.1)
+        else:
+            plt.suptitle(f"B = {B_fields[0]} mT")
+            plt.subplot(121)
+            plt.cla()
+            plt.plot(frequencies / u.MHz, np.reshape(R, len(frequencies)))
+            plt.xlabel("Qubit pulse frequency [MHz]")
+            plt.ylabel(r"$R=\sqrt{I^2 + Q^2}$ [V]")
+            plt.subplot(122)
+            plt.cla()
+            plt.plot(frequencies / u.MHz, np.reshape(phase, len(frequencies)))
+            plt.xlabel("Qubit pulse frequency [MHz]")
+            plt.ylabel("Phase [rad]")
+            plt.tight_layout()
+            plt.pause(0.1)
diff --git a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/12_randomized_benchmarking_single_qubit.py b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/12_randomized_benchmarking_single_qubit.py
new file mode 100644
index 000000000..a206925ea
--- /dev/null
+++ b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/12_randomized_benchmarking_single_qubit.py
@@ -0,0 +1,508 @@
+"""
+        SINGLE QUBIT RANDOMIZED BENCHMARKING (for gates >= 40ns)
+The program consists in playing random sequences of Clifford gates and measuring the state of the resonator afterwards.
+Each random sequence is derived on the FPGA for the maximum depth (specified as an input) and played for each depth
+asked by the user (the sequence is truncated to the desired depth). Each truncated sequence ends with the recovery gate,
+found at each step thanks to a preloaded lookup table (Cayley table), that will bring the qubit back to its ground state.
+
+If the readout has been calibrated and is good enough, then state discrimination can be applied to only return the state
+of the qubit. Otherwise, the 'I' and 'Q' quadratures are returned.
+Each sequence is played n_avg times for averaging. A second averaging is performed by playing different random sequences.
+
+The data is then post-processed to extract the single-qubit gate fidelity and error per gate
+.
+Prerequisites:
+    - Having found the resonance frequency of the resonator coupled to the qubit under study (resonator_spectroscopy).
+    - Having calibrated qubit pi pulse (x180) by running qubit, spectroscopy, rabi_chevron, power_rabi and updated the config.
+    - Having the qubit frequency perfectly calibrated (ramsey).
+    - (optional) Having calibrated the readout (readout_frequency, amplitude, duration_optimization IQ_blobs) for better SNR.
+"""
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool
+from qualang_tools.plot import interrupt_on_close
+from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
+from qualang_tools.voltage_gates import VoltageGateSequence
+from qualang_tools.addons.variables import assign_variables_to_element
+from macros import RF_reflectometry_macro, DC_current_sensing_macro
+from scipy.optimize import curve_fit
+import matplotlib.pyplot as plt
+
+
+##############################
+# Program-specific variables #
+##############################
+# Number of of averages for each random sequence
+n_avg = 50
+num_of_sequences = 50  # Number of random sequences
+max_circuit_depth = 1000  # Maximum circuit depth
+delta_clifford = 10  #  Play each sequence with a depth step equals to 'delta_clifford - Must be > 0
+assert (max_circuit_depth / delta_clifford).is_integer(), "max_circuit_depth / delta_clifford must be an integer."
+
+seed = 34553  # Pseudo-random number generator seed
+
+# Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)
+state_discrimination = False
+
+ge_threshold = 0.155  # arbitrary atm, in V
+B_field = 0  # Predetermined magnetic field value if not set externally
+
+
+# seq = VoltageGateSequence(config, ["P1_sticky", "P2_sticky"])
+seq = OPX_virtual_gate_sequence(config, ["P1_sticky", "P2_sticky"])
+seq.add_points("initialization", level_init, duration_init)
+# Idle is when RB sequence takes place, duration is overridden with calculated sequence timing
+seq.add_points("idle", level_manip, duration_manip)
+seq.add_points("readout", level_readout, duration_readout)
+
+# Time in ns for RB sequence to execute play_sequence(), will be buffer after ramping to idle
+# Note, with low max_circuit_depth, the delay before readout will increase slightly
+RB_delay = 92
+
+###################################
+# Helper functions and QUA macros #
+###################################
+
+# List of recovery gates from the lookup table
+inv_gates = [int(np.where(c1_table[i, :] == 0)[0][0]) for i in range(24)]
+
+
+def power_law(power, a, b, p):
+    return a * (p**power) + b
+
+
+def generate_sequence():
+    cayley = declare(int, value=c1_table.flatten().tolist())
+    inv_list = declare(int, value=inv_gates)
+    current_state = declare(int)
+    step = declare(int)
+    sequence = declare(int, size=max_circuit_depth + 1)
+    inv_gate = declare(int, size=max_circuit_depth + 1)
+    i = declare(int)
+    rand = Random(seed=seed)
+
+    assign(current_state, 0)
+    with for_(i, 0, i < max_circuit_depth, i + 1):
+        assign(step, rand.rand_int(24))
+        assign(current_state, cayley[current_state * 24 + step])
+        assign(sequence[i], step)
+        assign(inv_gate[i], inv_list[current_state])
+
+    return sequence, inv_gate
+
+
+def play_sequence(sequence_list, depth):
+    i = declare(int)
+    with for_(i, 0, i <= depth, i + 1):
+        with switch_(sequence_list[i], unsafe=True):
+            with case_(0):
+                wait(x180_len // 4, "qubit")
+            with case_(1):
+                play("x180", "qubit")
+            with case_(2):
+                play("y180", "qubit")
+            with case_(3):
+                play("y180", "qubit")
+                play("x180", "qubit")
+            with case_(4):
+                play("x90", "qubit")
+                play("y90", "qubit")
+            with case_(5):
+                play("x90", "qubit")
+                play("-y90", "qubit")
+            with case_(6):
+                play("-x90", "qubit")
+                play("y90", "qubit")
+            with case_(7):
+                play("-x90", "qubit")
+                play("-y90", "qubit")
+            with case_(8):
+                play("y90", "qubit")
+                play("x90", "qubit")
+            with case_(9):
+                play("y90", "qubit")
+                play("-x90", "qubit")
+            with case_(10):
+                play("-y90", "qubit")
+                play("x90", "qubit")
+            with case_(11):
+                play("-y90", "qubit")
+                play("-x90", "qubit")
+            with case_(12):
+                play("x90", "qubit")
+            with case_(13):
+                play("-x90", "qubit")
+            with case_(14):
+                play("y90", "qubit")
+            with case_(15):
+                play("-y90", "qubit")
+            with case_(16):
+                play("-x90", "qubit")
+                play("y90", "qubit")
+                play("x90", "qubit")
+            with case_(17):
+                play("-x90", "qubit")
+                play("-y90", "qubit")
+                play("x90", "qubit")
+            with case_(18):
+                play("x180", "qubit")
+                play("y90", "qubit")
+            with case_(19):
+                play("x180", "qubit")
+                play("-y90", "qubit")
+            with case_(20):
+                play("y180", "qubit")
+                play("x90", "qubit")
+            with case_(21):
+                play("y180", "qubit")
+                play("-x90", "qubit")
+            with case_(22):
+                play("x90", "qubit")
+                play("y90", "qubit")
+                play("x90", "qubit")
+            with case_(23):
+                play("-x90", "qubit")
+                play("y90", "qubit")
+                play("-x90", "qubit")
+
+
+# Macro to calculate exact duration of generated sequence at a given depth
+def generate_sequence_time(sequence_list, depth):
+    j = declare(int)
+    duration = declare(int)
+    assign(duration, 0)  # Ensures duration is reset to 0 for every depth calculated
+    with for_(j, 0, j <= depth, j + 1):
+        with switch_(sequence_list[j], unsafe=True):
+            with case_(0):
+                # wait(x180_len // 4, "qubit")
+                assign(duration, duration + x180_len)
+            with case_(1):
+                # play("x180", "qubit")
+                assign(duration, duration + x180_len)
+            with case_(2):
+                # play("y180", "qubit")
+                assign(duration, duration + y180_len)
+            with case_(3):
+                # play("y180", "qubit")
+                # play("x180", "qubit")
+                assign(duration, duration + y180_len + x180_len)
+            with case_(4):
+                # play("x90", "qubit")
+                # play("y90", "qubit")
+                assign(duration, duration + x90_len + y90_len)
+            with case_(5):
+                # play("x90", "qubit")
+                # play("-y90", "qubit")
+                assign(duration, duration + x90_len + minus_y90_len)
+            with case_(6):
+                # play("-x90", "qubit")
+                # play("y90", "qubit")
+                assign(duration, duration + minus_x90_len + y90_len)
+            with case_(7):
+                # play("-x90", "qubit")
+                # play("-y90", "qubit")
+                assign(duration, duration + minus_x90_len + minus_y90_len)
+            with case_(8):
+                # play("y90", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + y90_len + x90_len)
+            with case_(9):
+                # play("y90", "qubit")
+                # play("-x90", "qubit")
+                assign(duration, duration + y90_len + minus_x90_len)
+            with case_(10):
+                # play("-y90", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + minus_y90_len + x90_len)
+            with case_(11):
+                # play("-y90", "qubit")
+                # play("-x90", "qubit")
+                assign(duration, duration + minus_y90_len + minus_x90_len)
+            with case_(12):
+                # play("x90", "qubit")
+                assign(duration, duration + x90_len)
+            with case_(13):
+                # play("-x90", "qubit")
+                assign(duration, duration + minus_x90_len)
+            with case_(14):
+                # play("y90", "qubit")
+                assign(duration, duration + y90_len)
+            with case_(15):
+                # play("-y90", "qubit")
+                assign(duration, duration + minus_y90_len)
+            with case_(16):
+                # play("-x90", "qubit")
+                # play("y90", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + minus_x90_len + y90_len + x90_len)
+            with case_(17):
+                # play("-x90", "qubit")
+                # play("-y90", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + minus_x90_len + minus_y90_len + x90_len)
+            with case_(18):
+                # play("x180", "qubit")
+                # play("y90", "qubit")
+                assign(duration, duration + x180_len + y90_len)
+            with case_(19):
+                # play("x180", "qubit")
+                # play("-y90", "qubit")
+                assign(duration, duration + x180_len + minus_y90_len)
+            with case_(20):
+                # play("y180", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + y180_len + x90_len)
+            with case_(21):
+                # play("y180", "qubit")
+                # play("-x90", "qubit")
+                assign(duration, duration + y180_len + minus_x90_len)
+            with case_(22):
+                # play("x90", "qubit")
+                # play("y90", "qubit")
+                # play("x90", "qubit")
+                assign(duration, duration + x90_len + y90_len + x90_len)
+            with case_(23):
+                # play("-x90", "qubit")
+                # play("y90", "qubit")
+                # play("-x90", "qubit")
+                assign(duration, duration + minus_x90_len + y90_len + minus_x90_len)
+    return duration
+
+
+###################
+# The QUA program #
+###################
+with program() as rb:
+    depth = declare(int)  # QUA variable for the varying depth
+    depth_target = declare(int)  # QUA variable for the current depth (changes in steps of delta_clifford)
+    # QUA variable to store the last Clifford gate of the current sequence which is replaced by the recovery gate
+    saved_gate = declare(int)
+    m = declare(int)  # QUA variable for the loop over random sequences
+    n = declare(int)  # QUA variable for the averaging loop
+    I = declare(fixed)  # QUA variable for the 'I' quadrature
+    Q = declare(fixed)  # QUA variable for the 'Q' quadrature
+    state = declare(bool)  # QUA variable for state discrimination
+    sequence_time = declare(int)  # QUA variable for RB sequence duration for a given depth
+    dc_signal = declare(fixed)  # QUA variable for the measured dc signal
+    # Ensure that the result variables are assigned to the measurement elements
+    assign_variables_to_element("tank_circuit", I, Q)
+    assign_variables_to_element("TIA", dc_signal)
+
+    # The relevant streams
+    m_st = declare_stream()
+    if state_discrimination:
+        state_st = declare_stream()
+    else:
+        I_st = declare_stream()
+        Q_st = declare_stream()
+
+    with for_(m, 0, m < num_of_sequences, m + 1):  # QUA for_ loop over the random sequences
+
+        sequence_list, inv_gate_list = generate_sequence()  # Generate the random sequence of length max_circuit_depth
+
+        assign(depth_target, 1)  # Initialize the current depth to 1
+
+        with for_(depth, 1, depth <= max_circuit_depth, depth + 1):  # Loop over the depths
+
+            # Replacing the last gate in the sequence with the sequence's inverse gate
+            # The original gate is saved in 'saved_gate' and is being restored at the end
+            assign(saved_gate, sequence_list[depth])
+            assign(sequence_list[depth], inv_gate_list[depth - 1])
+
+            # Only played the depth corresponding to target_depth
+            with if_(depth == depth_target):
+
+                # Assign sequence_time to duration of idle step for generated sequence "m" at a given depth
+                assign(sequence_time, generate_sequence_time(sequence_list, depth))
+
+                with for_(n, 0, n < n_avg, n + 1):  # Averaging loop
+
+                    # Define voltage steps
+                    seq.add_step(voltage_point_name="initialization")
+                    seq.add_step(
+                        voltage_point_name="idle", duration=sequence_time + RB_delay
+                    )  # Includes processing time for RB sequence
+                    seq.add_step(voltage_point_name="readout", duration=readout_len)
+                    seq.add_compensation_pulse(duration=duration_compensation_pulse)
+
+                    wait((duration_init * u.ns), "qubit")
+                    play_sequence(sequence_list, depth)
+
+                    wait(
+                        (duration_init * u.ns) + (sequence_time >> 2) + (RB_delay >> 2),
+                        "tank_circuit",
+                        "TIA",
+                    )  # Includes calculated RB duration as well as RB sequence processing time
+                    I, Q, I_st, Q_st = RF_reflectometry_macro(I=I, Q=Q)
+                    dc_signal, dc_signal_st = DC_current_sensing_macro(dc_signal=dc_signal)
+
+                    # Make sure you updated the ge_threshold and angle if you want to use state discrimination
+                    # Save the results to their respective streams
+                    if state_discrimination:
+                        save(state, state_st)
+                    else:
+                        save(I, I_st)
+                        save(Q, Q_st)
+
+                    seq.ramp_to_zero()
+
+                # Go to the next depth
+                assign(depth_target, depth_target + delta_clifford)
+            # Reset the last gate of the sequence back to the original Clifford gate
+            # (that was replaced by the recovery gate at the beginning)
+            assign(sequence_list[depth], saved_gate)
+        # Save the counter for the progress bar
+        save(m, m_st)
+
+    with stream_processing():
+        m_st.save("iteration")
+        if state_discrimination:
+            # saves a 2D array of depth and random pulse sequences in order to get error bars along the random sequences
+            state_st.boolean_to_int().buffer(n_avg).map(FUNCTIONS.average()).buffer(
+                max_circuit_depth / delta_clifford
+            ).buffer(num_of_sequences).save("state")
+            # returns a 1D array of averaged random pulse sequences vs depth of circuit for live plotting
+            state_st.boolean_to_int().buffer(n_avg).map(FUNCTIONS.average()).buffer(
+                max_circuit_depth / delta_clifford
+            ).average().save("state_avg")
+        else:
+            I_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(max_circuit_depth / delta_clifford).buffer(
+                num_of_sequences
+            ).save("I")
+            Q_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(max_circuit_depth / delta_clifford).buffer(
+                num_of_sequences
+            ).save("Q")
+            I_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(max_circuit_depth / delta_clifford).average().save(
+                "I_avg"
+            )
+            Q_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(max_circuit_depth / delta_clifford).average().save(
+                "Q_avg"
+            )
+
+            # DC current sensing
+            dc_signal_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(max_circuit_depth / delta_clifford).buffer(
+                num_of_sequences
+            ).save("dc_signal")
+            dc_signal_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(
+                max_circuit_depth / delta_clifford
+            ).average().save("dc_signal_avg")
+
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+simulate = True
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=40_000)  # In clock cycles = 4ns
+    job = qmm.simulate(config, rb, simulation_config)
+
+    plt.figure()
+    job.get_simulated_samples().con1.plot()
+    # Get the waveform report
+    samples = job.get_simulated_samples()
+    waveform_report = job.get_simulated_waveform_report()
+    waveform_report.create_plot(samples, plot=True, save_path=None)
+    plt.legend("")
+    plt.show()
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(rb)
+    # Get results from QUA program
+    if state_discrimination:
+        results = fetching_tool(job, data_list=["state_avg", "iteration"], mode="live")
+    else:
+        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"], mode="live")
+    # Live plotting
+    fig = plt.figure()
+    interrupt_on_close(fig, job)  # Interrupts the job when closing the figure
+    # data analysis
+    x = np.arange(1, max_circuit_depth + 0.1, delta_clifford)
+    while results.is_processing():
+        # data analysis
+        if state_discrimination:
+            state_avg, iteration = results.fetch_all()
+            value_avg = state_avg
+        else:
+            I, Q, iteration = results.fetch_all()
+            value_avg = I
+
+        print(job.execution_report())
+        # Progress bar
+        progress_counter(iteration, num_of_sequences, start_time=results.get_start_time())
+        # Plot averaged values
+        plt.cla()
+        plt.plot(x, value_avg, marker=".")
+        plt.xlabel("Number of Clifford gates")
+        plt.ylabel("Sequence Fidelity")
+        plt.title("Single qubit RB")
+        plt.pause(0.1)
+
+    # At the end of the program, fetch the non-averaged results to get the error-bars
+    if state_discrimination:
+        results = fetching_tool(job, data_list=["state"])
+        state = results.fetch_all()[0]
+        value_avg = np.mean(state, axis=0)
+        error_avg = np.std(state, axis=0)
+    else:
+        results = fetching_tool(job, data_list=["I", "Q"])
+        I, Q = results.fetch_all()
+        value_avg = np.mean(I, axis=0)
+        error_avg = np.std(I, axis=0)
+
+    # data analysis
+    pars, cov = curve_fit(
+        f=power_law,
+        xdata=x,
+        ydata=value_avg,
+        p0=[0.5, 0.5, 0.9],
+        bounds=(-np.inf, np.inf),
+        maxfev=2000,
+    )
+    stdevs = np.sqrt(np.diag(cov))
+
+    print("#########################")
+    print("### Fitted Parameters ###")
+    print("#########################")
+    print(f"A = {pars[0]:.3} ({stdevs[0]:.1}), B = {pars[1]:.3} ({stdevs[1]:.1}), p = {pars[2]:.3} ({stdevs[2]:.1})")
+    print("Covariance Matrix")
+    print(cov)
+
+    one_minus_p = 1 - pars[2]
+    r_c = one_minus_p * (1 - 1 / 2**1)
+    r_g = r_c / 1.875  # 1.875 is the average number of gates in clifford operation
+    r_c_std = stdevs[2] * (1 - 1 / 2**1)
+    r_g_std = r_c_std / 1.875
+
+    print("#########################")
+    print("### Useful Parameters ###")
+    print("#########################")
+    print(
+        f"Error rate: 1-p = {np.format_float_scientific(one_minus_p, precision=2)} ({stdevs[2]:.1})\n"
+        f"Clifford set infidelity: r_c = {np.format_float_scientific(r_c, precision=2)} ({r_c_std:.1})\n"
+        f"Gate infidelity: r_g = {np.format_float_scientific(r_g, precision=2)}  ({r_g_std:.1})"
+    )
+
+    # Plots
+    plt.figure()
+    plt.errorbar(x, value_avg, yerr=error_avg, marker=".")
+    plt.plot(x, power_law(x, *pars), linestyle="--", linewidth=2)
+    plt.xlabel("Number of Clifford gates")
+    plt.ylabel("Sequence Fidelity")
+    plt.title("Single qubit RB")
+
+    # np.savez("rb_values", value)
+    # Close the quantum machines at the end in order to put all flux biases to 0 so that the fridge doesn't heat-up
+    qm.close()
diff --git a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/README.md b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/README.md
index c687fd560..fe698a946 100644
--- a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/README.md
+++ b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/README.md
@@ -39,15 +39,21 @@ is not saturated, correct for DC offsets.
 7. **Pauli Spin Blockade search** - Apply a triangle scan through the fast line of the bias-tess and on top of the charge stability map acquisition for finding the PSB readout point.
     * [Using the QDAC2 triggered by the OPX](07_PSB_search_qdac2_triggered.py)
     * [Using another external DC source](07_PSB_search_external_dc_source.py)
-8. [Qubit spectroscopy](08_qubit_spectroscopy.py) - Apply a strong and long qubit pulse and sweep its frequency over a wide range to identify the resonance. The magnetic field can also be swept to acquire the EDSR line and measure the g-factor.
+8. **Qubit Spectroscopy**
+    [Qubit spectroscopy](08a_qubit_spectroscopy.py) - Apply a strong and long qubit pulse and sweep its frequency over a wide range to identify the resonance. The magnetic field can also be swept to acquire the EDSR line and measure the g-factor.
+    **Chirp** -Alternative to qubit spectroscopy
+    [Chirp](8b_qubit_spectroscopy_with_chirp.py) -Allows user to define chirp duration and rate to sweep the IF and LO frequencies for quick determination of resonance.
 9. **Rabi chevron** - Measure the Rabi chevron by sweeping the qubit pulse frequency and duration. 
     * [Using real-time QUA](09a_rabi_chevron_qua.py) - Allows to sweep the pulse duration from 16ns and in steps of at least 4ns. There is no limit in the maximum pulse length or the number of points in the sweep.
     * [Using the baking tool](09b_rabi_chevron_baking.py) - Allows to sweep the pulse duration from 0ns and in steps of at least 1ns. Since the pulses must be loaded beforehand (like for an AWG), there is a limit in the number of samples that the OPX can memorize (65k per pulse processor).
     * [Using a combination of real-time QUA and baking](09c_rabi_chevron_baking+qua.py) - Combine the previous two methods in order to perform long scans with 1ns resolution.
-11. [T1](10_T1.py) - Measures T1.
-12. **Ramsey chevron** - Perform a 2D sweep (detuning versus idle time) to acquire the Ramsey chevron pattern.
+10. [T1](10_T1.py) - Measures T1.
+11. **Ramsey chevron** - Perform a 2D sweep (detuning versus idle time) to acquire the Ramsey chevron pattern.
     * [Using real-time QUA](11a_ramsey_chevron_4ns.py) - Allows to sweep the pulse duration from 16ns and in steps of at least 4ns. There is no limit in the maximum pulse length or the number of points in the sweep.
     * [Using the baking tool](11b_ramsey_chevron_full_baking.py) - Bake the full sequence (pi/2 - idle - pi/2) to allow 1ns resolution for the pi/2 pulses and exchange interaction time.
+12. **Single Qubit Randomized Benchmarking** - Perform an RB measurement to determine single qubit gate fidelites.
+    *[Using real-time QUA](12_randomized_benchmarking_single_qubit.py) - Allows the user to define and run different depths 
+    of a randomly generated sequence of single qubit XY clifford gates, followed by the appropriate inverse gate. Maximum number of gates = 7000.
 
 ## Use Cases
 
diff --git a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/configuration.py b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/configuration.py
index c358ef865..a2ecba58c 100644
--- a/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/configuration.py
+++ b/Quantum-Control-Applications/Quantum-Dots/Single_Spin_EDSR/configuration.py
@@ -311,7 +311,7 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
 duration_init = 2500
 duration_manip = 1000
 duration_readout = readout_len + 100
-duration_compensation_pulse = 4 * u.us
+duration_compensation_pulse = 4 * u.us  # Note, may need to be increased when running long RB sequences
 
 # Step parameters
 step_length = 16  # in ns
@@ -344,6 +344,32 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
 cw_amp = 0.3  # in V
 cw_len = 100  # in ns
 
+# Chirp Pulse
+chirp_duration = 1000  # in clock cycles
+chirp_rate = 5000
+chirp_units = "Hz/nsec"
+chirp_amp = 0.05  # 0.3
+processing_time = 196  # time in ns for chirp to be calculated
+
+# RB Gate specifics
+x180_len = 100
+x180_amp = 0.1  # 0.35
+
+x90_len = x180_len
+x90_amp = x180_amp / 2
+
+minus_x90_len = x180_len
+minus_x90_amp = -x90_amp
+
+y180_len = x180_len
+y180_amp = x180_amp
+
+y90_len = x180_len
+y90_amp = y180_amp / 2
+
+minus_y90_len = y180_len
+minus_y90_amp = -y90_amp
+
 #############################################
 #                  Config                   #
 #############################################
@@ -459,6 +485,13 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
                 "pi": "pi_pulse",
                 "pi_half": "pi_half_pulse",
                 "gauss": "gaussian_pulse",
+                "chirp": "chirp_pulse",
+                "x180": "x180_pulse",
+                "x90": "x90_pulse",
+                "-x90": "minus_x90_pulse",
+                "y180": "y180_pulse",
+                "y90": "y90_pulse",
+                "-y90": "minus_y90_pulse",
             },
         },
         "tank_circuit": {
@@ -526,6 +559,14 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
                 "Q": "zero_wf",
             },
         },
+        "chirp_pulse": {
+            "operation": "control",
+            "length": chirp_duration,
+            "waveforms": {
+                "I": "chirp_wf",
+                "Q": "zero_wf",
+            },
+        },
         "gaussian_pulse": {
             "operation": "control",
             "length": gaussian_length,
@@ -573,6 +614,54 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
             },
             "digital_marker": "ON",
         },
+        "x180_pulse": {
+            "operation": "control",
+            "length": x180_len,
+            "waveforms": {
+                "I": "x180_I_wf",
+                "Q": "x180_Q_wf",
+            },
+        },
+        "x90_pulse": {
+            "operation": "control",
+            "length": x90_len,
+            "waveforms": {
+                "I": "x90_I_wf",
+                "Q": "x90_Q_wf",
+            },
+        },
+        "minus_x90_pulse": {
+            "operation": "control",
+            "length": x90_len,
+            "waveforms": {
+                "I": "-x90_I_wf",
+                "Q": "-x90_Q_wf",
+            },
+        },
+        "y180_pulse": {
+            "operation": "control",
+            "length": y180_len,
+            "waveforms": {
+                "I": "y180_I_wf",
+                "Q": "y180_Q_wf",
+            },
+        },
+        "y90_pulse": {
+            "operation": "control",
+            "length": y90_len,
+            "waveforms": {
+                "I": "y90_I_wf",
+                "Q": "y90_Q_wf",
+            },
+        },
+        "minus_y90_pulse": {
+            "operation": "control",
+            "length": y90_len,
+            "waveforms": {
+                "I": "-y90_I_wf",
+                "Q": "-y90_Q_wf",
+            },
+        },
     },
     "waveforms": {
         "P1_step_wf": {"type": "constant", "sample": P1_step_amp},
@@ -588,6 +677,20 @@ def add_points(self, name: str, coordinates: list, duration: int) -> None:
         "reflect_wf": {"type": "constant", "sample": reflectometry_readout_amp},
         "const_wf": {"type": "constant", "sample": cw_amp},
         "zero_wf": {"type": "constant", "sample": 0.0},
+        "chirp_wf": {"type": "constant", "sample": chirp_amp},
+        "zero_wf": {"type": "constant", "sample": 0.0},
+        "x180_I_wf": {"type": "constant", "sample": x180_amp},
+        "x180_Q_wf": {"type": "constant", "sample": 0.0},
+        "x90_I_wf": {"type": "constant", "sample": x90_amp},
+        "x90_Q_wf": {"type": "constant", "sample": 0.0},
+        "-x90_I_wf": {"type": "constant", "sample": minus_x90_amp},
+        "-x90_Q_wf": {"type": "constant", "sample": 0.0},
+        "y180_I_wf": {"type": "constant", "sample": 0.0},
+        "y180_Q_wf": {"type": "constant", "sample": y180_amp},
+        "y90_I_wf": {"type": "constant", "sample": 0.0},
+        "y90_Q_wf": {"type": "constant", "sample": y90_amp},
+        "-y90_I_wf": {"type": "constant", "sample": 0.0},
+        "-y90_Q_wf": {"type": "constant", "sample": minus_y90_amp},
     },
     "digital_waveforms": {
         "ON": {"samples": [(1, 0)]},
diff --git a/Quantum-Control-Applications/README.md b/Quantum-Control-Applications/README.md
index 796e2a115..59e99c129 100644
--- a/Quantum-Control-Applications/README.md
+++ b/Quantum-Control-Applications/README.md
@@ -12,6 +12,7 @@ These files showcase various experiments that can be done on a several flux-tuna
 * [SWAP spectroscopy improved with predistortion digital filters](https://github.com/qua-platform/qua-libs/tree/main/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Transmons/Use%20Case%201%20-%20Two%20qubit%20gate%20optimization%20with%20cryoscope#two-qubit-swap-spectroscopy-improved-with-pre-distortion-digital-filters).
 * [Two-qubit randomized benchmarking](https://github.com/qua-platform/qua-libs/tree/main/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Transmons/Use%20Case%202%20-%20Two-Qubit-Randomized-Benchmarking#two-qubit-randomized-benchmarking).
 * [Two-qubit cross-entropy benchmarking](https://github.com/qua-platform/qua-libs/tree/main/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Transmons/Use%20Case%203%20-%20Two-Qubit%20Cross-Entropy%20Benchmarking).
+* [Single- and Two-Qubit State and Process Tomography](https://github.com/qua-platform/qua-libs/tree/main/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use%20Case%204%20-%20Single-%20and%20Two-Qubit%20State%20and%20Process%20Tomography)
 
 ### [Single fixed frequency transmon](https://github.com/qua-platform/qua-libs/tree/main/Quantum-Control-Applications/Superconducting/Single-Fixed-Transmon#single-fixed-transmon-superconducting-qubit)
 These files showcase various experiments that can be done on a single fixed-frequency transmon.
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01a_manual_mixer_calibration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01a_manual_mixer_calibration.py
new file mode 100644
index 000000000..71b3b9162
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01a_manual_mixer_calibration.py
@@ -0,0 +1,114 @@
+"""
+        MIXER CALIBRATION
+The program is designed to play a continuous single tone to calibrate an IQ mixer. To do this, connect the mixer's
+output to a spectrum analyzer. Adjustments for the DC offsets, gain, and phase must be made manually.
+
+If you have access to the API for retrieving data from the spectrum analyzer, you can utilize the commented lines below
+to semi-automate the process.
+
+Before proceeding to the next node, take the following steps:
+    - Update the DC offsets in the configuration at: config/controllers/"con1"/analog_outputs.
+    - Modify the DC gain and phase for the IQ signals in the configuration, under either:
+      mixer_qubit_g & mixer_qubit_g or mixer_resonator_g & mixer_resonator_g.
+"""
+
+from qm import QuantumMachinesManager
+from qm.qua import *
+from configuration import *
+
+###################
+# The QUA program #
+###################
+element = "rr1"
+
+with program() as cw_output:
+    with infinite_loop_():
+        # It is best to calibrate LO leakage first and without any power played (cf. note below)
+        play("cw" * amp(0), element)
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+qm = qmm.open_qm(config)
+
+job = qm.execute(cw_output)
+
+# When done, the halt command can be called and the offsets can be written directly into the config file.
+
+# job.halt()
+
+# These are the 2 commands used to correct for mixer imperfections. The first is used to set the DC of the `I` and `Q`
+# channels to compensate for the LO leakage. The 2nd command is used to correct for the phase and amplitude mismatches
+# between the channels.
+# The output of the IQ Mixer should be connected to a spectrum analyzer and values should be chosen as to minimize the
+# unwanted peaks.
+# If python can read the output of the spectrum analyzer, then this process can be automated and the correct values can
+# be found using an optimization method such as Nelder-Mead:
+# https://docs.scipy.org/doc/scipy/reference/optimize.minimize-neldermead.html
+
+# qm.set_output_dc_offset_by_element('rr1', ('I', 'Q'), (-0.001, 0.003))
+# qm.set_mixer_correction('mixer_resonator', int(resonator_IF_q1), int(resonator_LO), IQ_imbalance(0.015, 0.01))
+
+# Note that the LO leakage (DC Offset) depends on the 'I' & 'Q' powers, it is advised to run this step with no input power.
+# This will ensure that there is no LO leakage while the pulses are not played in the case where the is no switch.
+# This can be achieved by changing the line above to `play("cw" * amp(0), "qubit")`
+
+# Automatic LO leakage correction
+# centers = [0.5, 0]
+# span = 0.1
+#
+# fig1 = plt.figure()
+# for n in range(3):
+#     offset_i = np.linspace(centers[0] - span, centers[0] + span, 21)
+#     offset_q = np.linspace(centers[1] - span, centers[1] + span, 31)
+#     lo_leakage = np.zeros((len(offset_q), len(offset_i)))
+#     for i in range(len(offset_i)):
+#         for q in range(len(offset_q)):
+#             qm.set_output_dc_offset_by_element(element, ("I", "Q"), (offset_i[i], offset_q[q]))
+#             sleep(0.01)
+#             # Write functions to extract the lo leakage from the spectrum analyzer
+#             # lo_leakage[q][i] =
+#     minimum = np.argwhere(lo_leakage == np.min(lo_leakage))[0]
+#     centers = [offset_i[minimum[0]], offset_q[minimum[1]]]
+#     span = span / 10
+#     plt.subplot(131)
+#     plt.pcolor(offset_i, offset_q, lo_leakage.transpose())
+#     plt.xlabel("I offset [V]")
+#     plt.ylabel("Q offset [V]")
+#     plt.title(f"Minimum at (I={centers[0]:.3f}, Q={centers[1]:.3f}) = {lo_leakage[minimum[0]][minimum[1]]:.1f} dBm")
+# plt.suptitle(f"LO leakage correction for {element}")
+#
+# print(f"For {element}, I offset is {centers[0]} and Q offset is {centers[1]}")
+#
+# # Automatic image cancellation
+# centers = [0.5, 0]
+# span = [0.2, 0.5]
+#
+# fig2 = plt.figure()
+# for n in range(3):
+#     gain = np.linspace(centers[0] - span, centers[0] + span, 21)
+#     phase = np.linspace(centers[1] - span, centers[1] + span, 31)
+#     image = np.zeros((len(phase), len(gain)))
+#     for g in range(len(gain)):
+#         for p in range(len(phase)):
+#             qm.set_mixer_correction(
+#                 config["elements"][element]["mixInputs"]["mixer"],
+#                 int(config["elements"][element]["intermediate_frequency"]),
+#                 int(config["elements"][element]["mixInputs"]["lo_frequency"]),
+#                 IQ_imbalance(gain[g], phase[p]),
+#             )
+#             sleep(0.01)
+#             # Write functions to extract the image from the spectrum analyzer
+#             # image[q][i] =
+#     minimum = np.argwhere(image == np.min(image))[0]
+#     centers = [gain[minimum[0]], phase[minimum[1]]]
+#     span = (np.array(span) / 10).tolist()
+#     plt.subplot(131)
+#     plt.pcolor(gain, phase, image.transpose())
+#     plt.xlabel("Gain")
+#     plt.ylabel("Phase imbalance [rad]")
+#     plt.title(f"Minimum at (I={centers[0]:.3f}, Q={centers[1]:.3f}) = {image[minimum[0]][minimum[1]]:.1f} dBm")
+# plt.suptitle(f"Image cancellation for {element}")
+#
+# print(f"For {element}, gain is {centers[0]} and phase is {centers[1]}")
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01b_octave_clock_and_calibration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01b_octave_clock_and_calibration.py
new file mode 100644
index 000000000..f95ac1f8e
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01b_octave_clock_and_calibration.py
@@ -0,0 +1,32 @@
+"""
+This file is used to configure the Octave's clock and do the automatic calibration.
+"""
+
+from qm import QuantumMachinesManager
+from qm.octave import ClockMode
+from configuration import *
+
+
+# Configure the Octave according to the elements settings and calibrate
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, octave=octave_config, log_level="ERROR")
+qm = qmm.open_qm(config)
+
+##################
+# Clock settings #
+##################
+qm.octave.set_clock("octave1", clock_mode=ClockMode.Internal)
+# If using external LO change this line to one of the following:
+#     qm.octave.set_clock("octave1", clock_mode=ClockMode.External_10MHz)
+#     qm.octave.set_clock("octave1", clock_mode=ClockMode.External_100MHz)
+#     qm.octave.set_clock("octave1", clock_mode=ClockMode.External_1000MHz)
+
+##################
+# Calibration #
+##################
+calibration = True
+
+if calibration:
+    elements = ["rr1", "rr2", "q1_xy", "q2_xy"]
+    for element in elements:
+        print("-" * 37 + f" Calibrates {element}")
+        qm.calibrate_element(element)
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01_time_of_flight.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/02_time_of_flight.py
similarity index 100%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/01_time_of_flight.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/02_time_of_flight.py
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/02_resonator_spectroscopy_single.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_single.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/02_resonator_spectroscopy_single.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_single.py
index d2a86a9cf..af0709002 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/02_resonator_spectroscopy_single.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_single.py
@@ -17,7 +17,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager, SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 from qualang_tools.results import fetching_tool
 from qualang_tools.loops import from_array
 import matplotlib.pyplot as plt
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_multiplexed.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_multiplexed.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_multiplexed.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_multiplexed.py
index 89728eb8a..40d50cf40 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/03_resonator_spectroscopy_multiplexed.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_multiplexed.py
@@ -17,7 +17,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager, SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 from qualang_tools.results import progress_counter, fetching_tool
 from qualang_tools.plot import interrupt_on_close
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_vs_amplitude.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_resonator_spectroscopy_vs_amplitude.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_vs_amplitude.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_resonator_spectroscopy_vs_amplitude.py
index 34a9fe028..9309b9213 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/04_resonator_spectroscopy_vs_amplitude.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_resonator_spectroscopy_vs_amplitude.py
@@ -21,7 +21,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager, SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 from qualang_tools.results import progress_counter, fetching_tool
 from qualang_tools.plot import interrupt_on_close
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_qubit_spectroscopy.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_qubit_spectroscopy.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_qubit_spectroscopy.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_qubit_spectroscopy.py
index 1f26deaa8..62f83b7fa 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/05_qubit_spectroscopy.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_qubit_spectroscopy.py
@@ -26,7 +26,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager, SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 from qualang_tools.results import progress_counter, fetching_tool
 from qualang_tools.plot import interrupt_on_close
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_rabi_chevron.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_rabi_chevron.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_rabi_chevron.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_rabi_chevron.py
index 1d898e3e4..a1a9e5327 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/06_rabi_chevron.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_rabi_chevron.py
@@ -18,7 +18,7 @@
 from qm import QuantumMachinesManager
 from qm.qua import *
 from qm import SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_power_rabi.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_power_rabi.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_power_rabi.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_power_rabi.py
index f3c6e3316..5fb30a6c7 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/07_power_rabi.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_power_rabi.py
@@ -18,7 +18,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_ramsey_chevron.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09_ramsey_chevron.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_ramsey_chevron.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09_ramsey_chevron.py
index 3241c4404..806ffdaf5 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/08_ramsey_chevron.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09_ramsey_chevron.py
@@ -16,7 +16,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09a_readout_optimization_freq.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10a_readout_optimization_freq.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09a_readout_optimization_freq.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10a_readout_optimization_freq.py
index b66a6f5bc..31f666b34 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09a_readout_optimization_freq.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10a_readout_optimization_freq.py
@@ -18,7 +18,7 @@
 from qm import QuantumMachinesManager
 from qm.qua import *
 from qm import SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09b_readout_optimization_amp.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10b_readout_optimization_amp.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09b_readout_optimization_amp.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10b_readout_optimization_amp.py
index 6a357de12..5d22b7ad0 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09b_readout_optimization_amp.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10b_readout_optimization_amp.py
@@ -17,7 +17,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09c_readout_optimization_duration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10c_readout_optimization_duration.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09c_readout_optimization_duration.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10c_readout_optimization_duration.py
index 11147e7fd..5a4d9ab03 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09c_readout_optimization_duration.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10c_readout_optimization_duration.py
@@ -4,7 +4,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.results import fetching_tool, progress_counter
 import math
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09d_readout_weight_optimization.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10d_readout_weight_optimization.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09d_readout_weight_optimization.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10d_readout_weight_optimization.py
index 2553ab598..50528ec57 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/09d_readout_weight_optimization.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10d_readout_weight_optimization.py
@@ -4,7 +4,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.results import fetching_tool, progress_counter
 from qualang_tools.plot import interrupt_on_close
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10_IQ_blobs.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_IQ_blobs.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10_IQ_blobs.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_IQ_blobs.py
index 586a2e32c..d41e5a2ba 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/10_IQ_blobs.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_IQ_blobs.py
@@ -19,7 +19,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.results import fetching_tool
 from macros import qua_declaration, multiplexed_readout, active_reset
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_T1.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T1.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_T1.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T1.py
index fa5c9b98d..7d32e5f1c 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/11_T1.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T1.py
@@ -16,7 +16,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.results import fetching_tool
 from qualang_tools.plot import interrupt_on_close
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T2echo.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2echo.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T2echo.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2echo.py
index 9228e7200..047222849 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/12_T2echo.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2echo.py
@@ -4,7 +4,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.results import fetching_tool
 from qualang_tools.plot import interrupt_on_close
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2ramsey.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14_T2ramsey.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2ramsey.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14_T2ramsey.py
index ab57feedb..494638197 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/13_T2ramsey.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14_T2ramsey.py
@@ -16,7 +16,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14a_ac_stark_calibration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15a_ac_stark_calibration.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14a_ac_stark_calibration.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15a_ac_stark_calibration.py
index a8b10a362..f72626873 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14a_ac_stark_calibration.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15a_ac_stark_calibration.py
@@ -22,7 +22,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14b_drag_calibration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15b_drag_calibration.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14b_drag_calibration.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15b_drag_calibration.py
index 5959fb2f8..b8c204019 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/14b_drag_calibration.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15b_drag_calibration.py
@@ -18,7 +18,7 @@
 
 from qm import QuantumMachinesManager, SimulationConfig
 from qm.qua import *
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qualang_tools.loops import from_array
 from qualang_tools.results import fetching_tool, progress_counter
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15_allxy.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16_allxy.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15_allxy.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16_allxy.py
index bb3f1284f..04f129043 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/15_allxy.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16_allxy.py
@@ -4,7 +4,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.results import fetching_tool
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16a_single_qubit_RB.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_single_qubit_RB.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16a_single_qubit_RB.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_single_qubit_RB.py
index ba5636840..a1e59033b 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16a_single_qubit_RB.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_single_qubit_RB.py
@@ -5,7 +5,7 @@
 from qm.qua import *
 from qm import QuantumMachinesManager
 from qm import SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 import numpy as np
 from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16b_single_qubit_RB_interleaved.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_single_qubit_RB_interleaved.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16b_single_qubit_RB_interleaved.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_single_qubit_RB_interleaved.py
index d768056a5..735c25811 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/16b_single_qubit_RB_interleaved.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_single_qubit_RB_interleaved.py
@@ -5,7 +5,7 @@
 from qm.qua import *
 from qm import QuantumMachinesManager
 from qm import SimulationConfig
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 import numpy as np
 from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_CR_time_rabi_1q_QST.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18a_CR_time_rabi_1q_QST.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_CR_time_rabi_1q_QST.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18a_CR_time_rabi_1q_QST.py
index ba934e30d..c58200340 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17a_CR_time_rabi_1q_QST.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18a_CR_time_rabi_1q_QST.py
@@ -18,7 +18,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_cancelCR_time_rabi_1q_QST.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_cancelCR_time_rabi_1q_QST.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_cancelCR_time_rabi_1q_QST.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_cancelCR_time_rabi_1q_QST.py
index 12bc8b79f..d4b6933b3 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17b_cancelCR_time_rabi_1q_QST.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_cancelCR_time_rabi_1q_QST.py
@@ -18,7 +18,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17c_echoCR_time_rabi_1q_QST.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_echoCR_time_rabi_1q_QST.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17c_echoCR_time_rabi_1q_QST.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_echoCR_time_rabi_1q_QST.py
index 7026ef5d1..b64abd132 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/17c_echoCR_time_rabi_1q_QST.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_echoCR_time_rabi_1q_QST.py
@@ -18,7 +18,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18a_CR_calib_unit_hamiltonian_tomography.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19a_CR_calib_unit_hamiltonian_tomography.py
similarity index 100%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18a_CR_calib_unit_hamiltonian_tomography.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19a_CR_calib_unit_hamiltonian_tomography.py
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_CR_calib_cr_drive_amplitude.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19b_CR_calib_cr_drive_amplitude.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_CR_calib_cr_drive_amplitude.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19b_CR_calib_cr_drive_amplitude.py
index 54151d64e..37d198adb 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18b_CR_calib_cr_drive_amplitude.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19b_CR_calib_cr_drive_amplitude.py
@@ -39,7 +39,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_CR_calib_cr_drive_phase.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19c_CR_calib_cr_drive_phase.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_CR_calib_cr_drive_phase.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19c_CR_calib_cr_drive_phase.py
index 10a967f9f..1b0d58bc6 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18c_CR_calib_cr_drive_phase.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19c_CR_calib_cr_drive_phase.py
@@ -39,7 +39,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18d_CR_calib_cr_cancel_phase.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19d_CR_calib_cr_cancel_phase.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18d_CR_calib_cr_cancel_phase.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19d_CR_calib_cr_cancel_phase.py
index d3b8096b8..c0b679f8d 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18d_CR_calib_cr_cancel_phase.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19d_CR_calib_cr_cancel_phase.py
@@ -42,7 +42,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18e_CR_calib_cr_cancel_amplitude.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19e_CR_calib_cr_cancel_amplitude.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18e_CR_calib_cr_cancel_amplitude.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19e_CR_calib_cr_cancel_amplitude.py
index d7fbc4077..a77c99619 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18e_CR_calib_cr_cancel_amplitude.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19e_CR_calib_cr_cancel_amplitude.py
@@ -42,7 +42,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18f_CR_calib_cr_driven_ramsey_RCVersion.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19f_CR_calib_cr_driven_ramsey_RCVersion.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18f_CR_calib_cr_driven_ramsey_RCVersion.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19f_CR_calib_cr_driven_ramsey_RCVersion.py
index e8eb4ef93..0d14d2437 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/18f_CR_calib_cr_driven_ramsey_RCVersion.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19f_CR_calib_cr_driven_ramsey_RCVersion.py
@@ -37,7 +37,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19_CNOT.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/20_CNOT.py
similarity index 99%
rename from Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19_CNOT.py
rename to Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/20_CNOT.py
index c9d9a6cbf..b32780316 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/19_CNOT.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/20_CNOT.py
@@ -14,7 +14,7 @@
 
 from qm.qua import *
 from qm import QuantumMachinesManager
-from configuration_mw_fem import *
+from configuration import *
 import matplotlib.pyplot as plt
 from qm import SimulationConfig
 from qualang_tools.loops import from_array
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration.py
index 694f0256b..6f403422a 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration.py
@@ -33,8 +33,8 @@ def IQ_imbalance(g, phi):
 ######################
 # Network parameters #
 ######################
-qop_ip = "172.16.33.101"  # Write the QM router IP address
-cluster_name = "Cluster_83"  # Write your cluster_name if version >= QOP220
+qop_ip = "127.0.0.1"  # Write the QM router IP address
+cluster_name = None  # Write your cluster_name if version >= QOP220
 qop_port = None  # Write the QOP port if version < QOP220
 octave_config = None
 
@@ -419,8 +419,8 @@ def IQ_imbalance(g, phi):
         },
         "cr_cancel_c1t2": {
             "mixInputs": {
-                "I": ("con1", 1),
-                "Q": ("con1", 2),
+                "I": ("con1", 3),
+                "Q": ("con1", 4),
                 "lo_frequency": qubit_LO_q1,
                 "mixer": "mixer_qubit_q1",
             },
@@ -433,8 +433,8 @@ def IQ_imbalance(g, phi):
         },
         "cr_cancel_c2t1": {
             "mixInputs": {
-                "I": ("con1", 3),
-                "Q": ("con1", 4),
+                "I": ("con1", 1),
+                "Q": ("con1", 2),
                 "lo_frequency": qubit_LO_q2,
                 "mixer": "mixer_qubit_q2",
             },
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem.py
index 3162663f8..8bee15c16 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem.py
@@ -32,8 +32,8 @@ def IQ_imbalance(g, phi):
 ######################
 # Network parameters #
 ######################
-qop_ip = "172.16.33.116"  # Write the QM router IP address
-cluster_name = "Beta_8"  # Write your cluster_name if version >= QOP220
+qop_ip = "127.0.0.1"  # Write the QM router IP address
+cluster_name = None  # Write your cluster_name if version >= QOP220
 qop_port = None  # Write the QOP port if version < QOP220
 octave_config = None
 
@@ -483,8 +483,8 @@ def IQ_imbalance(g, phi):
         },
         "cr_cancel_c1t2": {
             "mixInputs": {
-                "I": (con, fem, 1),
-                "Q": (con, fem, 2),
+                "I": (con, fem, 3),
+                "Q": (con, fem, 4),
                 "lo_frequency": qubit_LO_q1,
                 "mixer": "mixer_qubit_q1",
             },
@@ -497,8 +497,8 @@ def IQ_imbalance(g, phi):
         },
         "cr_cancel_c2t1": {
             "mixInputs": {
-                "I": (con, fem, 3),
-                "Q": (con, fem, 4),
+                "I": (con, fem, 1),
+                "Q": (con, fem, 2),
                 "lo_frequency": qubit_LO_q2,
                 "mixer": "mixer_qubit_q2",
             },
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem_and_octave.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem_and_octave.py
index 2aa7f813b..ec13bd540 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem_and_octave.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_lf_fem_and_octave.py
@@ -33,8 +33,8 @@ def IQ_imbalance(g, phi):
 ######################
 # Network parameters #
 ######################
-qop_ip = "172.16.33.116"  # Write the QM router IP address
-cluster_name = "Beta_8"  # Write your cluster_name if version >= QOP220
+qop_ip = "127.0.0.1"  # Write the QM router IP address
+cluster_name = None  # Write your cluster_name if version >= QOP220
 qop_port = None  # Write the QOP port if version < QOP220
 
 ############################
@@ -467,7 +467,7 @@ def IQ_imbalance(g, phi):
             },
         },
         "cr_cancel_c1t2": {
-            "RF_inputs": {"port": ("octave1", 2)},
+            "RF_inputs": {"port": ("octave1", 3)},
             "intermediate_frequency": cr_cancel_IF_c1t2,  # in Hz
             "operations": {
                 "cw": "const_pulse",
@@ -476,7 +476,7 @@ def IQ_imbalance(g, phi):
             },
         },
         "cr_cancel_c2t1": {
-            "RF_inputs": {"port": ("octave1", 3)},
+            "RF_inputs": {"port": ("octave1", 2)},
             "intermediate_frequency": cr_cancel_IF_c2t1,  # in Hz
             "operations": {
                 "cw": "const_pulse",
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_mw_fem.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_mw_fem.py
index d6650015f..a31279e57 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_mw_fem.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_mw_fem.py
@@ -18,8 +18,8 @@
 ######################
 # Network parameters #
 ######################
-qop_ip = "172.16.33.116"  # Write the QM router IP address
-cluster_name = "Beta_8"  # Write your cluster_name if version >= QOP220
+qop_ip = "127.0.0.1"  # Write the QM router IP address
+cluster_name = None  # Write your cluster_name if version >= QOP220
 qop_port = None  # Write the QOP port if version < QOP220
 octave_config = None
 
diff --git a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_with_octave.py b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_with_octave.py
index 346470270..99e357bde 100644
--- a/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_with_octave.py
+++ b/Quantum-Control-Applications/Superconducting/Two-Fixed-Coupled-Transmons/configuration_with_octave.py
@@ -33,8 +33,8 @@ def IQ_imbalance(g, phi):
 ######################
 # Network parameters #
 ######################
-qop_ip = "172.16.33.101"  # Write the QM router IP address
-cluster_name = "Cluster_81"  # Write your cluster_name if version >= QOP220
+qop_ip = "127.0.0.1"  # Write the QM router IP address
+cluster_name = None  # Write your cluster_name if version >= QOP220
 qop_port = None  # Write the QOP port if version < QOP220
 
 ############################
@@ -399,7 +399,7 @@ def IQ_imbalance(g, phi):
             },
         },
         "cr_cancel_c1t2": {
-            "RF_inputs": {"port": ("octave1", 2)},
+            "RF_inputs": {"port": ("octave1", 3)},
             "intermediate_frequency": cr_cancel_IF_c1t2,  # in Hz
             "operations": {
                 "cw": "const_pulse",
@@ -408,7 +408,7 @@ def IQ_imbalance(g, phi):
             },
         },
         "cr_cancel_c2t1": {
-            "RF_inputs": {"port": ("octave1", 3)},
+            "RF_inputs": {"port": ("octave1", 2)},
             "intermediate_frequency": cr_cancel_IF_c2t1,  # in Hz
             "operations": {
                 "cw": "const_pulse",
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/LICENSE.txt b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/LICENSE.txt
new file mode 100644
index 000000000..70ca978ec
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/LICENSE.txt	
@@ -0,0 +1,11 @@
+                              Copyright 2024. Fermi Research Alliance, LLC.
+					All rights reserved.
+
+Copyright 2024. Fermi Research Alliance, LLC. This software was produced under U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Superconducting Quantum Materials and Systems (SQMS) Center, under the contract No. DE-AC02-07CH11359, for Fermi National Accelerator Laboratory (Fermilab), which is operated by Fermi Research Alliance, LLC for the U.S. Department of Energy. The U.S. Government has rights to use, reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR FERMI RESEARCH ALLIANCE, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.  If software is modified to produce derivative works, such modified software should be clearly marked, so as not to confuse it with the version available from Fermilab.
+ 
+Additionally, redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+• Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+• Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+• Neither the name of Fermi Research Alliance, LLC, Fermi National Accelerator Laboratory, Fermilab, the U.S. Government, nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY FERMI RESEARCH ALLIANCE, LLC AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL FERMI RESEARCH ALLIANCE, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
\ No newline at end of file
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/PMatrix2.pkl b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/PMatrix2.pkl
new file mode 100644
index 000000000..aa40bc0f4
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/PMatrix2.pkl differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/README.md b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/README.md
new file mode 100644
index 000000000..1edc47c4c
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/README.md	
@@ -0,0 +1,185 @@
+# Single- and Two-Qubit State and Process Tomography
+
+_Authors and experimenters: Nicholas Bornman, Paul Heidler and Ziwen Huang_
+
+*Important note: the code in this folder is the code that was used for running these measurements on a pair of flux-tunable*
+*qubits coupled with a flux-tunable coupler. The code may need to be modified for a different piece of hardware*
+*and software environment.*
+
+This use-case has four potential applications, some of which may not be applicable/necessary, namely
+1. Single-qubit state tomography,
+2. Single-qubit process tomography,
+3. Two-qubit state tomography, and
+4. Two-qubit process tomography
+
+The most complex applications are obviously the two-qubit cases, which is what we will focus on below.
+The single-qubit cases follow trivially from them.
+
+## 1. Experimental set-up
+
+<img align="right" src="setup.png" alt="setup" width="400"/>
+
+
+### 1.1 Experimental set-up
+The chip consists of 9 2D flux-tunable transmons arranged into a 3x3 grid, with nearest neighbours connected via flux tunable lines.
+
+For the purposes of this experiment, we focus on two neighbouring flux-tunable transmons (tuned using SQUID non-linear inductors and nearby flux lines) coupled to individual readout resonators with a common feedline, as well as the tunable coupler between them. This system is isolated from the rest of the chip by parking neighbouring qubits far away via external DC sources. See the figure on the right.
+
+The qubits are driven with analog IQ signal output pairs generated by the OPX+ and up-converted with an Octave. Their individual fluxes, as well as the flux of the tunable coupler, are tuned using three additional analog outputs from the OPX+.
+
+The readout resonators are driven with an analog IF signals pair generated by the OPX+ and up-converted with an Octave.
+
+The transmitted signal is measured by the OPX after down-conversion with the Octave, and the response of each resonator extracted by employing a multiplexed readout scheme.
+
+Note: single qubit tomography could be carried out using this chip by focusing on one of the two qubits and parking the second qubit far away.
+
+
+### 1.2 Pre-requisite calibrations
+
+
+For each resonator and qubit:
+
+* The first calibration is to find the readout resonator frequency and readout pulse amplitude using 2D resonator spectroscopy for each resonator.
+* With the frequencies and readout amplitudes set, perform spectroscopy on each resonator while sweeping the flux bias of the corresponding coupled qubit. The resulting 2D "cosine" plots show the resonators' frequencies as functions of the flux biases. From this, the fluxes at the qubits' sweet spots, at which point the qubits are most insensitive to flux noise, can be set. We want to set the flux lines to the qubits' flux insensitive points by default.
+* Next, we calibrate the qubits' $\pi/2$ and $\pi$ pulses. This can generally be done by first performing qubit spectroscopy to get a rough idea of the qubits' frequencies. Next, simple Rabi chevron measurements can be carried out, where both the $pi$ pulse amplitudes and frequencies are scanned while the flux bias is set to the qubit flux insensitive point found previously. The qubits' frequencies can be fine-tuned using Ramsey measurements, and the $\pi$ pulses calibrated further by performing Rabi error amplification experiments.
+* (optional) If the qubit pulses are relatively short (< 100ns), then the use of optimal DRAG waveforms can enhance the fidelity of the gates.
+* A two-state readout discriminator, one for each qubit, needs to be calibrated by performing readout fidelity "IQ blob" measurements, and the discriminator thresholds as well as the discriminator rotation angles, set.
+* Find the DC flux offset value of the tunable coupler such that the coupling between the two qubits in question, is off. We want the coupling between the qubits to be off, by default.
+
+Note that the state preparation and/or measurement "parts" of the state/process tomography, in the current tomography schemes, do not make use of the flux tunability of the qubits nor that of the tunable coupler, beyond parking the qubits at constant DC offsets. This is due to the fact that only separable two-qubit operations (i.e. two separate $\pi$ or $\pi/2$ pulses) are applied to prepare and/or measure the qubits' state. However, the user may wish to analyse an entangled state/pulse 'gate' operation'. In such cases, the required qua commands, which may use different constant DC fluxes or flux pulses as part of the subroutine to create the entanglement, need to be entered into the "prepare_state"/"analysed_process" functions of the scripts.
+
+
+## 2. [The configuration](configuration.py)
+
+
+### 2.1 The elements
+The configuration consists of 4 elements:
+* `rr1` and `rr2` send readout pulses and measure the transmitted signal for each qubit
+* `q1_xy` and `q2_xy` send drive XY pulses to each qubit
+* `q1_z` and `q2_z` send control Z pulses and DC voltage to the flux line of each qubit
+* `tc12` sends control pulses and DC voltage to the flux loop of the tunable coupler
+
+
+### 2.2 The operations, pulse and waveforms
+
+The readout resonators `rr1` and `rr2` each have a 'readout' operation, each of which is a constant pulse with calibrated readout lengths and amplitudes.
+
+The qubits elements `q1_xy` and `q2_xy` have several operations defined that correspond to X/2, Y/2, -X/2, -Y/2, X and Y gates. The envelopes are Gaussian (although this can be changed by the user), with optional Gaussian DRAG corrections, parameters (amplitudes, lengths, Gaussian sigmas, etc.), that need to be calibrated.
+
+The qubits' flux elements, `q1_z` and `q2_z`, have a single "zero_flux" operation that is not used. Instead, the constant DC offsets (which are configured in the controllers section) which set the qubits at their flux sweet spots throughout an entire QUA program, need to be calibrated and the values assigned to 'max_frequency_point_q1' and 'max_frequency_point_q2'.
+
+The tunable coupler `tc12`, similarly, has a single "zero_flux" operation that is not used, but rather the constant DC offset value for the corresponding flux line needs to be calibrated such that the coupling between qubits 1 and 2 is turned off, and the value of this offset assigned is to 't12_coupling_off_flux_value'.
+
+
+
+## 3. The QUA program(s) and data processing
+
+The QUA programs we will discuss below, fall into two distinct parts:
+1. Two-qubit state tomography, and
+2. Two-qubit process tomography
+
+Note, as mentioned, that we have single matching single qubit tomography scripts; with an understanding of the two-qubit case established, the single qubit case should be relatively trivial.
+
+For detailed explanations of tomography and detailed accompanying simulations, the interested reader is encouraged to consult [[1]](#1). Only the salient points are discussed below.
+
+
+### 3.1 [Two-qubit state tomography](two-qubit-state-tomography.py)
+
+
+The idea of this script is to perform standard state tomography on a two-qubit state chosen and prepared by the user, by performing a tomographically-complete set of two-qubit measurements on an ensemble of identically-prepared such states.
+
+We assume that the two qubits are always initialised in their ground states. The users needs to determine which sequence of QUA play statements are needed in order to prepare the desired two-qubit state in a noiseless, ideal scenario, and modify the `prepare_state(qubit1, qubit2)` function accordingly. For example, starting with $\left | 0 \right \rangle \left | 0 \right \rangle$, a Y gate applied to the first (left) qubit and a -X/2 gate applied to the second qubit, should ideally produce the state $\frac{1}{\sqrt{2}} \left | 1 \right \rangle \left ( \left | 0 \right \rangle + i\left | 0 \right \rangle \right)$. Here, the body of the `prepare_state(qubit1, qubit2)` should be modified such that:
+
+```
+def prepare_state(qubit1, qubit2):
+    # write whatever QUA code you need in order to create the
+    # state to perform tomography on, from an initial
+    # ground state. For example, to create the |1> 1/sqrt(2)(|0>+i|1>)
+    # state
+    play("y180", f"q{qubit1}_xy")
+    play("-x90", f"q{qubit2}_xy")
+```
+
+Note that if one wishes to prepare a more complex, possible entangled two-qubit state, the configuration.py potentially needs to be modified to include additional elements and pulses. These pulses would need to be played within the `prepare_state` method.
+
+The user should also modify the `ideal_state` Python variable to a density matrix representation of the ideal two-qubit state (in order to compute the state fidelity). In the above example, this would be:
+
+```
+# (|1> 1/sqrt(2)(|0>+i|1>)) (<1| 1/sqrt(2)(<0|-i<1|))
+# where |0> = [1, 0]^T, |1> = [0, 1]^T
+ideal_state = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1 / 2, 1j / 2], [0, 0, -1j / 2, 1 / 2]])
+```
+
+This is all that should be required.
+
+When run, the QUA program prepares the state being examined and then performs one of 15 possible combinations of single qubit rotations (which corresponds to projective measurements in a number of different measurement bases). The calibrated discriminators are used to perform single shot measurements on the two qubit state, which corresponds to which of four possible computational two-qubit states were measured.
+
+This process is repeated for all 15 projective measurement bases, with the result streamed to a stream processing tag called "probs". The 15 possible Stokes parameters are computed from various combinations of the `probs` Python variable (see [1](#1) for details), which serve as the coefficients in the Pauli operator basis decomposition of any two-qubit density matrix.
+
+The final result is `rho`, a numpy array representing the experimental density matrix of the two-qubit state.
+
+
+### 3.2 [Two-qubit process tomography](two-qubit-process-tomography.py)
+
+
+In standard two-qubit process tomography, one is generally interested in studying the dynamics of an open quantum system described by a completely-positive linear map:
+
+$$ \epsilon(\rho) = \sum_{m,n=0}^{4^2 - 1} \chi_{mn} E_m \rho E_n^{\dagger}, $$
+
+where the $\chi$ process matrix completely characterises the mapping, and $\{ E_i \}$ are any fixed Hermitian basis for the space of linear operators acting on the two-qubit Hilbert space, and are usually taken to be the Pauli matrices ([[2]](#2) for details). In other words, one is interested in determining the superoperator $\chi$, from which one can potentially determine the unitary dynamics (in the ideal, noiseless case) and/or glean various information about potential noise channels that exist in the system. It should be noted that a common use case is when the dynamics are well represented by a two-qubit unitary operator $U$ acting on the system's Hilbert space, then $\chi$ assumes a special form such that the above reduces to $\epsilon(\rho) = U \rho U^{\dagger}$.
+
+Suppose that the process under study would in fact be such a two-qubit unitary operator $U$ in an ideal, noiseless scenario. Next, assume we had access to calibrated gates with which to carry out $U$ (which could themselves entail well-calibrated pulse sequences to effect a non-trivial entangling gate between the two qubits). In this case, the two-qubit process tomography script here prepares a tomographically-complete set of input states, plays the process pulse sequence to effect $U$, and then performs standard two-qubit state tomography on the result. The results are then post-processed to produce $\chi$.
+
+The user needs to determine the QUA statements subroutine that is needed to produce the unitary $U$ in an ideal scenario, and modify the `analysed_process(qubit1, qubit2)` function accordingly. For example, suppose that we wish to analyse the process of applying an X gate on qubit 1 and simultaneously a -Y/2 gate on qubit 2, such that the process is ideally $\epsilon(\rho) = U \rho U^{\dagger}$ where $U = X \otimes (-Y/2)$. The body of the `analysed_process(qubit1, qubit2)` should be modified as follows:
+
+```
+# subroutine to prepare desired gate/process
+def analysed_process(qubit1, qubit2):
+    # write whatever QUA code you need, here, in order to perform the
+    # desired process which we want to subject to tomography.
+    # For example, to analyse the X gate on qubit1 and the -Y90
+    # gate on qubit2:
+    play("x180", f"q{qubit1}_xy")
+    play("-y90", f"q{qubit2}_xy")
+```
+
+Note that if one wishes to prepare a more complex entanlging gate, the configuration.py would potentially need to be modified to include additional elements and/or pulses. The sequence of QUA instructions which would be needed in order to carry out the entangling operation would need to be coded within the `analysed_process` method.
+
+Optionally, the user could also modify the `ideal_gate` Python variable to a unitary matrix representation of the ideal two-qubit gate (in order to compute some fidelity metrics). By consulting the [helper funtions](https://github.com/bornman-nick/qua-libs/blob/main/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Transmons/Use%20Case%204%20-%20Single-%20and%20Two-Qubit%20State%20and%20Process%20Tomography/helper_functions.py) file, note that `func_F2(1)` is an X gate, and `func_F2(3)` a -Y/2 gate. With the above $U$ in mind, the `ideal_gate` is:
+
+```
+# X \otimes (-Y/2)
+ideal_gate = np.kron(func_F2(1), func_F2(3))
+```
+
+Note also that the absolute file path to the PMatrix2.pkl file should be assigned to the `file_location` variable, in order to load the `PMatrix` variable, a constant matrix which is needed in the process of determining $\chi$. If not specified, `PMatrix` will be computed from scratch, but this is somewhat computationally time-consuming.
+
+These are all the modifications required from the user.
+
+When run, the QUA program prepares one of the required 36 separable two-qubit input states, runs the `analysed_process` method to create an almost certainly noisy version of $U$, and then performs one of the required 36 possible combinations of single qubit `play` sequences to change the measurement basis (see [[1]](#1) for details). The calibrated discriminators are used to perform single shot measurements on the two-qubit state, with the resulting boolean variable (called `state`) indicating whether both qubits were in their ground states or not.
+
+This process is repeated for all 35 input states and 35 measurement bases, with the result streamed to a stream processing tag called "probs". The Python variable `probs` is post processed by a number of helper functions to produce `measurement_vector` (corresponding to $\lambda$ in Eq. 5 of [[2]](#2)), which is finally used to solve the necessary matrix equation for the variable `chi_matrix` in the script, which corresponds to $\chi$. The final result, $\chi$, as mentioned, contains all information about the channel's dynamics which may be necessary for further analysis by the user, and is finally plotted.
+
+Both the noiseless, theoretical $\chi$ as well as the $\chi$ measured from the two qubits of our experimental hardware, are plotted below for an implementation of the above simple example process, namely an X gate on qubit 1 and simultaneously a -Y/2 gate on qubit 2:
+
+|             Ideal                 |            Real                  |
+|:---------------------------------:|:--------------------------------:|
+| ![ideal_XmYon2](XmYon2_ideal.png) | ![real_XmYon2](XmYon2_real.png)  |
+
+
+Finally, we repeat the two-qubit process tomography for an entangling iSWAP operation (the calibration and details of which are beyond the scope of this readme, but the intersted reader is directed to [[3]](#3) for details):
+
+|             Ideal              |              Real             |
+|:------------------------------:|:-----------------------------:|
+| ![ideal_swap](iswap_ideal.png) | ![real_swap](iswap_real.png)  |
+
+
+One could use the experimental $\chi$ matrix to further characterise the dynamics, such as computing the process fidelity, or potentially determining the Kraus noise operators of the channel.
+
+## References
+
+<a id="1">[1]</a> Bornman, N. Quantum state and process tomography (2024). Git repository: https://github.com/bornman-nick/quantum-state-and-process-tomography
+
+<a id="2">[2]</a> Mohseni, M. and Rezakhani, A. T. and Lidar, D. A. Quantum Process Tomography: Resource Analysis of Different Strategies (2008). https://arxiv.org/pdf/quant-ph/0702131
+
+<a id="3">[3]</a> Feng, B. et. al. Fluxonium: An Alternative Qubit Platform for High-Fidelity Operations (2022). https://link.aps.org/doi/10.1103/PhysRevLett.129.010502
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_ideal.png b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_ideal.png
new file mode 100644
index 000000000..2f1a7f520
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_ideal.png differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_real.png b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_real.png
new file mode 100644
index 000000000..d0992345f
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/XmYon2_real.png differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/configuration.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/configuration.py
new file mode 100644
index 000000000..42e109708
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/configuration.py	
@@ -0,0 +1,676 @@
+"""
+Octave configuration working for QOP222 and qm-qua==1.1.5 and newer.
+"""
+
+from pathlib import Path
+import numpy as np
+from set_octave import OctaveUnit, octave_declaration
+from qualang_tools.config.waveform_tools import drag_gaussian_pulse_waveforms
+from qualang_tools.units import unit
+
+
+#######################
+# AUXILIARY FUNCTIONS #
+#######################
+u = unit(coerce_to_integer=True)
+
+######################
+# Network parameters #
+######################
+qop_ip = "192.168.5.7"  # Write the QM router IP address
+cluster_name = "my_cluster"  # Write your cluster_name if version >= QOP220
+qop_port = None  # Write the QOP port if version < QOP220
+
+# Path to save data
+save_dir = Path().absolute() / "QM" / "INSTALLATION" / "data"
+
+############################
+# Set octave configuration #
+############################
+
+# The Octave port is 11xxx, where xxx are the last three digits of the Octave internal IP that can be accessed from
+# the OPX admin panel if you QOP version is >= QOP220. Otherwise, it is 50 for Octave1, then 51, 52 and so on.
+octave_1 = OctaveUnit("octave1", qop_ip, port=11253, con="con1")
+# octave_2 = OctaveUnit("octave2", qop_ip, port=11051, con="con1")
+
+# If the control PC or local network is connected to the internal network of the QM router (port 2 onwards)
+# or directly to the Octave (without QM the router), use the local octave IP and port 80.
+# octave_ip = "192.168.88.X"
+# octave_1 = OctaveUnit("octave1", octave_ip, port=80, con="con1")
+
+
+# Add the octaves
+octaves = [octave_1]
+# Configure the Octaves
+octave_config = octave_declaration(octaves)
+
+#####################
+# OPX configuration #
+#####################
+
+#############################################
+#                  Qubits                   #
+#############################################
+qubit_LO_q1 = 4.75 * u.GHz
+qubit_LO_q2 = qubit_LO_q1
+
+# Qubits IF
+qubit_IF_q1 = (-90) * u.MHz
+qubit_IF_q2 = (-158) * u.MHz
+
+# Relaxation time
+T1_q1 = int(50 * u.us)
+T1_q2 = int(50 * u.us)
+thermalization_time = 10 * max(T1_q1, T1_q2)
+
+# CW pulse parameter
+const_len = 50000
+const_amp = 200 * u.mV
+
+# Pi pulse parameters
+pi_len = 60 * u.ns
+pi_sigma = pi_len / 5
+pi_amp_q1 = 0.14
+pi_amp_q2 = 0.15
+
+# DRAG coefficients
+drag_coef_q1 = 1.6
+drag_coef_q2 = 1.6
+anharmonicity_q1 = (-180) * u.MHz
+anharmonicity_q2 = (-180) * u.MHz
+AC_stark_detuning_q1 = 0 * u.MHz
+AC_stark_detuning_q2 = 0 * u.MHz
+
+
+# DRAG waveforms
+
+
+# X pulses
+
+x180_wf_q1, x180_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q1, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1)
+)
+x180_I_wf_q1 = x180_wf_q1
+x180_Q_wf_q1 = x180_der_wf_q1
+
+x180_wf_q2, x180_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2)
+)
+x180_I_wf_q2 = x180_wf_q2
+x180_Q_wf_q2 = x180_der_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+# X/2 pulses
+
+x90_wf_q1, x90_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q1 / 2, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1)
+)
+x90_I_wf_q1 = x90_wf_q1
+x90_Q_wf_q1 = x90_der_wf_q1
+
+x90_wf_q2, x90_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q2 / 2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2)
+)
+x90_I_wf_q2 = x90_wf_q2
+x90_Q_wf_q2 = x90_der_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+# -X/2 pulses
+
+minus_x90_wf_q1, minus_x90_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(
+        -pi_amp_q1 / 2, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1
+    )
+)
+minus_x90_I_wf_q1 = minus_x90_wf_q1
+minus_x90_Q_wf_q1 = minus_x90_der_wf_q1
+
+minus_x90_wf_q2, minus_x90_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(
+        -pi_amp_q2 / 2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2
+    )
+)
+minus_x90_I_wf_q2 = minus_x90_wf_q2
+minus_x90_Q_wf_q2 = minus_x90_der_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+# Y pulses
+
+y180_wf_q1, y180_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q1, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1)
+)
+y180_I_wf_q1 = (-1) * y180_der_wf_q1
+y180_Q_wf_q1 = y180_wf_q1
+
+y180_wf_q2, y180_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2)
+)
+y180_I_wf_q2 = (-1) * y180_der_wf_q2
+y180_Q_wf_q2 = y180_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+# Y/2 pulses
+
+y90_wf_q1, y90_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q1 / 2, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1)
+)
+y90_I_wf_q1 = (-1) * y90_der_wf_q1
+y90_Q_wf_q1 = y90_wf_q1
+
+y90_wf_q2, y90_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(pi_amp_q2 / 2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2)
+)
+y90_I_wf_q2 = (-1) * y90_der_wf_q2
+y90_Q_wf_q2 = y90_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+# -Y/2 pulses
+
+minus_y90_wf_q1, minus_y90_der_wf_q1 = np.array(
+    drag_gaussian_pulse_waveforms(
+        -pi_amp_q1 / 2, pi_len, pi_sigma, drag_coef_q1, anharmonicity_q1, AC_stark_detuning_q1
+    )
+)
+minus_y90_I_wf_q1 = (-1) * minus_y90_der_wf_q1
+minus_y90_Q_wf_q1 = minus_y90_wf_q1
+
+minus_y90_wf_q2, minus_y90_der_wf_q2 = np.array(
+    drag_gaussian_pulse_waveforms(
+        -pi_amp_q2 / 2, pi_len, pi_sigma, drag_coef_q2, anharmonicity_q2, AC_stark_detuning_q2
+    )
+)
+minus_y90_I_wf_q2 = (-1) * minus_y90_der_wf_q2
+minus_y90_Q_wf_q2 = minus_y90_wf_q2
+# No DRAG when alpha=0, it's just a gaussian.
+
+
+##########################################
+#               Flux line                #
+##########################################
+flux_settle_time = 40 * u.ns
+
+const_flux_len = 3 * u.us
+
+max_frequency_point_q1 = 0.0
+max_frequency_point_q2 = 0.0
+
+min_frequency_point_q1 = 0.25
+min_frequency_point_q2 = 0.2490
+
+t12_coupling_off_flux_value = 0.1220
+t12_coupling_on_flux_value = 0
+
+# Resonator frequency versus flux fit parameters according to resonator_spec_vs_flux
+# amplitude * np.cos(2 * np.pi * frequency * x + phase) + offset (Hz, 1/V, degrees, Hz)
+amplitude_fit_q1, frequency_fit_q1, phase_fit_q1, offset_fit_q1 = [0, 0, 0, 0]
+amplitude_fit_q2, frequency_fit_q2, phase_fit_q2, offset_fit_q2 = [0, 0, 0, 0]
+
+
+#############################################
+#                Resonators                 #
+#############################################
+resonator_LO = 7.2 * u.GHz
+
+# Resonators IF
+resonator_IF_q1 = int(64 * u.MHz)
+resonator_IF_q2 = int(-111 * u.MHz)
+
+# Readout pulse parameters
+readout_len = 3000
+readout_amp_q1 = 0.01
+readout_amp_q2 = 0.01
+
+# TOF and depletion time
+time_of_flight = 240  # must be a multiple of 4
+depletion_time = 50 * u.us
+
+opt_weights = False
+if opt_weights:
+    from qualang_tools.config.integration_weights_tools import convert_integration_weights
+
+    weights_q1 = np.load("optimal_weights_q1.npz")
+    opt_weights_real_q1 = convert_integration_weights(weights_q1["weights_real"])
+    opt_weights_minus_imag_q1 = convert_integration_weights(weights_q1["weights_minus_imag"])
+    opt_weights_imag_q1 = convert_integration_weights(weights_q1["weights_imag"])
+    opt_weights_minus_real_q1 = convert_integration_weights(weights_q1["weights_minus_real"])
+
+    weights_q2 = np.load("optimal_weights_q2.npz")
+    opt_weights_real_q2 = convert_integration_weights(weights_q2["weights_real"])
+    opt_weights_minus_imag_q2 = convert_integration_weights(weights_q2["weights_minus_imag"])
+    opt_weights_imag_q2 = convert_integration_weights(weights_q2["weights_imag"])
+    opt_weights_minus_real_q2 = convert_integration_weights(weights_q2["weights_minus_real"])
+
+else:
+
+    opt_weights_real_q1 = [(1.0, readout_len)]
+    opt_weights_minus_imag_q1 = [(0.0, readout_len)]
+    opt_weights_imag_q1 = [(0.0, readout_len)]
+    opt_weights_minus_real_q1 = [(-1.0, readout_len)]
+
+    opt_weights_real_q2 = [(1.0, readout_len)]
+    opt_weights_minus_imag_q2 = [(0.0, readout_len)]
+    opt_weights_imag_q2 = [(0.0, readout_len)]
+    opt_weights_minus_real_q2 = [(-1.0, readout_len)]
+
+# state discrimination
+
+rotation_angle_q1 = ((0) / 180) * np.pi
+rotation_angle_q2 = ((0) / 180) * np.pi
+ge_threshold_q1 = 0.0
+ge_threshold_q2 = 0.0
+
+
+#############################################
+#                  Config                   #
+#############################################
+
+config = {
+    "version": 1,
+    "controllers": {
+        "con1": {
+            "analog_outputs": {
+                1: {"offset": 0.0},  # I readout line
+                2: {"offset": 0.0},  # Q readout line
+                3: {"offset": max_frequency_point_q1},  # qubit1 Z
+                4: {"offset": max_frequency_point_q2},  # qubit2 Z
+                5: {"offset": 0.0},  # I qubit1 XY
+                6: {"offset": 0.0},  # Q qubit1 XY
+                7: {"offset": 0.0},  # I qubit2 XY
+                8: {"offset": 0.0},  # Q qubit2 XY
+                9: {"offset": 0.0},  # N/A
+                10: {"offset": t12_coupling_off_flux_value},  # qubit1-qubit2 coupler
+            },
+            "digital_outputs": {
+                1: {},
+                3: {},
+                5: {},
+                7: {},
+                9: {},
+            },
+            "analog_inputs": {
+                1: {"offset": 0.0, "gain_db": 0},  # I from down-conversion
+                2: {"offset": 0.0, "gain_db": 0},  # Q from down-conversion
+            },
+        },
+    },
+    "elements": {
+        "rr1": {
+            "RF_inputs": {"port": ("octave1", 1)},
+            "RF_outputs": {"port": ("octave1", 1)},
+            "intermediate_frequency": resonator_IF_q1,
+            "operations": {
+                "readout": "readout_pulse_q1",
+            },
+            "digitalInputs": {
+                "switch": {
+                    "port": ("con1", 1),
+                    "delay": 57,
+                    "buffer": 18,
+                },
+            },
+            "time_of_flight": time_of_flight,
+            "smearing": 0,
+        },
+        "rr2": {
+            "RF_inputs": {"port": ("octave1", 1)},
+            "RF_outputs": {"port": ("octave1", 1)},
+            "intermediate_frequency": resonator_IF_q2,
+            "operations": {
+                "readout": "readout_pulse_q2",
+            },
+            "digitalInputs": {
+                "switch": {
+                    "port": ("con1", 1),
+                    "delay": 57,
+                    "buffer": 18,
+                },
+            },
+            "time_of_flight": time_of_flight,
+            "smearing": 0,
+        },
+        "q1_xy": {
+            "RF_inputs": {"port": ("octave1", 3)},
+            "intermediate_frequency": qubit_IF_q1,
+            "operations": {
+                "cw": "const_pulse",
+                "x180": "x180_pulse_q1",
+                "x90": "x90_pulse_q1",
+                "-x90": "-x90_pulse_q1",
+                "y90": "y90_pulse_q1",
+                "y180": "y180_pulse_q1",
+                "-y90": "-y90_pulse_q1",
+            },
+            "digitalInputs": {
+                "switch": {
+                    "port": ("con1", 5),
+                    "delay": 57,
+                    "buffer": 18,
+                },
+            },
+        },
+        "q2_xy": {
+            "RF_inputs": {"port": ("octave1", 4)},
+            "intermediate_frequency": qubit_IF_q2,
+            "operations": {
+                "cw": "const_pulse",
+                "x180": "x180_pulse_q2",
+                "x90": "x90_pulse_q2",
+                "-x90": "-x90_pulse_q2",
+                "y90": "y90_pulse_q2",
+                "y180": "y180_pulse_q2",
+                "-y90": "-y90_pulse_q2",
+            },
+            "digitalInputs": {
+                "switch": {
+                    "port": ("con1", 7),
+                    "delay": 57,
+                    "buffer": 18,
+                },
+            },
+        },
+        "q1_z": {
+            "singleInput": {
+                "port": ("con1", 3),
+            },
+            "operations": {"zero": "zero_flux"},
+        },
+        "q2_z": {
+            "singleInput": {
+                "port": ("con1", 4),
+            },
+            "operations": {"zero": "zero_flux"},
+        },
+        "tc12": {
+            "singleInput": {
+                "port": ("con1", 10),
+            },
+            "operations": {
+                "zero": "zero_flux",
+            },
+        },
+    },
+    "octaves": {
+        "octave1": {
+            "RF_outputs": {
+                1: {
+                    "LO_frequency": resonator_LO,
+                    "LO_source": "internal",
+                    "output_mode": "always_on",
+                    "gain": 0,
+                },
+                3: {
+                    "LO_frequency": qubit_LO_q1,
+                    "LO_source": "internal",
+                    "output_mode": "always_on",  # "triggered_reversed",
+                    "gain": 0,
+                },
+                4: {
+                    "LO_frequency": qubit_LO_q2,
+                    "LO_source": "internal",
+                    "output_mode": "always_on",  # "triggered_reversed",
+                    "gain": 0,
+                },
+            },
+            "RF_inputs": {
+                1: {
+                    "LO_frequency": resonator_LO,
+                    "LO_source": "internal",
+                },
+            },
+            "connectivity": "con1",
+        }
+    },
+    "pulses": {
+        "const_pulse": {
+            "operation": "control",
+            "length": const_len,
+            "waveforms": {
+                "I": "const_wf",
+                "Q": "zero_wf",
+            },
+        },
+        "readout_pulse_q1": {
+            "operation": "measurement",
+            "length": readout_len,
+            "waveforms": {
+                "I": "readout_wf_q1",
+                "Q": "zero_wf",
+            },
+            "integration_weights": {
+                "cos": "cosine_weights",
+                "sin": "sine_weights",
+                "minus_sin": "minus_sine_weights",
+                "rotated_cos": "rotated_cosine_weights_q1",
+                "rotated_sin": "rotated_sine_weights_q1",
+                "rotated_minus_sin": "rotated_minus_sine_weights_q1",
+                "opt_cos": "opt_cosine_weights_q1",
+                "opt_sin": "opt_sine_weights_q1",
+                "opt_minus_sin": "opt_minus_sine_weights_q1",
+            },
+            "digital_marker": "ON",
+        },
+        "readout_pulse_q2": {
+            "operation": "measurement",
+            "length": readout_len,
+            "waveforms": {
+                "I": "readout_wf_q2",
+                "Q": "zero_wf",
+            },
+            "integration_weights": {
+                "cos": "cosine_weights",
+                "sin": "sine_weights",
+                "minus_sin": "minus_sine_weights",
+                "rotated_cos": "rotated_cosine_weights_q2",
+                "rotated_sin": "rotated_sine_weights_q2",
+                "rotated_minus_sin": "rotated_minus_sine_weights_q2",
+                "opt_cos": "opt_cosine_weights_q2",
+                "opt_sin": "opt_sine_weights_q2",
+                "opt_minus_sin": "opt_minus_sine_weights_q2",
+            },
+            "digital_marker": "ON",
+        },
+        "zero_flux": {
+            "operation": "control",
+            "length": const_flux_len,
+            "waveforms": {
+                "single": "zero_wf",
+            },
+        },
+        "x90_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "x90_I_wf_q1",
+                "Q": "x90_Q_wf_q1",
+            },
+        },
+        "x180_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "x180_I_wf_q1",
+                "Q": "x180_Q_wf_q1",
+            },
+        },
+        "-x90_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "minus_x90_I_wf_q1",
+                "Q": "minus_x90_Q_wf_q1",
+            },
+        },
+        "y90_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "y90_I_wf_q1",
+                "Q": "y90_Q_wf_q1",
+            },
+        },
+        "y180_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "y180_I_wf_q1",
+                "Q": "y180_Q_wf_q1",
+            },
+        },
+        "-y90_pulse_q1": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "minus_y90_I_wf_q1",
+                "Q": "minus_y90_Q_wf_q1",
+            },
+        },
+        "x90_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "x90_I_wf_q2",
+                "Q": "x90_Q_wf_q2",
+            },
+        },
+        "x180_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "x180_I_wf_q2",
+                "Q": "x180_Q_wf_q2",
+            },
+        },
+        "-x90_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "minus_x90_I_wf_q2",
+                "Q": "minus_x90_Q_wf_q2",
+            },
+        },
+        "y90_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "y90_I_wf_q2",
+                "Q": "y90_Q_wf_q2",
+            },
+        },
+        "y180_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "y180_I_wf_q2",
+                "Q": "y180_Q_wf_q2",
+            },
+        },
+        "-y90_pulse_q2": {
+            "operation": "control",
+            "length": pi_len,
+            "waveforms": {
+                "I": "minus_y90_I_wf_q2",
+                "Q": "minus_y90_Q_wf_q2",
+            },
+        },
+    },
+    "waveforms": {
+        "const_wf": {"type": "constant", "sample": const_amp},
+        "zero_wf": {"type": "constant", "sample": 0.0},
+        "x90_I_wf_q1": {"type": "arbitrary", "samples": x90_I_wf_q1.tolist()},
+        "x90_Q_wf_q1": {"type": "arbitrary", "samples": x90_Q_wf_q1.tolist()},
+        "x180_I_wf_q1": {"type": "arbitrary", "samples": x180_I_wf_q1.tolist()},
+        "x180_Q_wf_q1": {"type": "arbitrary", "samples": x180_Q_wf_q1.tolist()},
+        "minus_x90_I_wf_q1": {"type": "arbitrary", "samples": minus_x90_I_wf_q1.tolist()},
+        "minus_x90_Q_wf_q1": {"type": "arbitrary", "samples": minus_x90_Q_wf_q1.tolist()},
+        "y90_I_wf_q1": {"type": "arbitrary", "samples": y90_I_wf_q1.tolist()},
+        "y90_Q_wf_q1": {"type": "arbitrary", "samples": y90_Q_wf_q1.tolist()},
+        "y180_I_wf_q1": {"type": "arbitrary", "samples": y180_I_wf_q1.tolist()},
+        "y180_Q_wf_q1": {"type": "arbitrary", "samples": y180_Q_wf_q1.tolist()},
+        "minus_y90_I_wf_q1": {"type": "arbitrary", "samples": minus_y90_I_wf_q1.tolist()},
+        "minus_y90_Q_wf_q1": {"type": "arbitrary", "samples": minus_y90_Q_wf_q1.tolist()},
+        "readout_wf_q1": {"type": "constant", "sample": readout_amp_q1},
+        "x90_I_wf_q2": {"type": "arbitrary", "samples": x90_I_wf_q2.tolist()},
+        "x90_Q_wf_q2": {"type": "arbitrary", "samples": x90_Q_wf_q2.tolist()},
+        "x180_I_wf_q2": {"type": "arbitrary", "samples": x180_I_wf_q2.tolist()},
+        "x180_Q_wf_q2": {"type": "arbitrary", "samples": x180_Q_wf_q2.tolist()},
+        "minus_x90_I_wf_q2": {"type": "arbitrary", "samples": minus_x90_I_wf_q2.tolist()},
+        "minus_x90_Q_wf_q2": {"type": "arbitrary", "samples": minus_x90_Q_wf_q2.tolist()},
+        "y90_I_wf_q2": {"type": "arbitrary", "samples": y90_I_wf_q2.tolist()},
+        "y90_Q_wf_q2": {"type": "arbitrary", "samples": y90_Q_wf_q2.tolist()},
+        "y180_I_wf_q2": {"type": "arbitrary", "samples": y180_I_wf_q2.tolist()},
+        "y180_Q_wf_q2": {"type": "arbitrary", "samples": y180_Q_wf_q2.tolist()},
+        "minus_y90_I_wf_q2": {"type": "arbitrary", "samples": minus_y90_I_wf_q2.tolist()},
+        "minus_y90_Q_wf_q2": {"type": "arbitrary", "samples": minus_y90_Q_wf_q2.tolist()},
+        "readout_wf_q2": {"type": "constant", "sample": readout_amp_q2},
+    },
+    "digital_waveforms": {
+        "ON": {"samples": [(1, 0)]},
+    },
+    "integration_weights": {
+        "cosine_weights": {
+            "cosine": [(1.0, readout_len)],
+            "sine": [(0.0, readout_len)],
+        },
+        "sine_weights": {
+            "cosine": [(0.0, readout_len)],
+            "sine": [(1.0, readout_len)],
+        },
+        "minus_sine_weights": {
+            "cosine": [(0.0, readout_len)],
+            "sine": [(-1.0, readout_len)],
+        },
+        "rotated_cosine_weights_q2": {
+            "cosine": [(np.cos(rotation_angle_q2), readout_len)],
+            "sine": [(np.sin(rotation_angle_q2), readout_len)],
+        },
+        "rotated_sine_weights_q2": {
+            "cosine": [(-np.sin(rotation_angle_q2), readout_len)],
+            "sine": [(np.cos(rotation_angle_q2), readout_len)],
+        },
+        "rotated_minus_sine_weights_q2": {
+            "cosine": [(np.sin(rotation_angle_q2), readout_len)],
+            "sine": [(-np.cos(rotation_angle_q2), readout_len)],
+        },
+        "rotated_cosine_weights_q1": {
+            "cosine": [(np.cos(rotation_angle_q1), readout_len)],
+            "sine": [(np.sin(rotation_angle_q1), readout_len)],
+        },
+        "rotated_sine_weights_q1": {
+            "cosine": [(-np.sin(rotation_angle_q1), readout_len)],
+            "sine": [(np.cos(rotation_angle_q1), readout_len)],
+        },
+        "rotated_minus_sine_weights_q1": {
+            "cosine": [(np.sin(rotation_angle_q1), readout_len)],
+            "sine": [(-np.cos(rotation_angle_q1), readout_len)],
+        },
+        "opt_cosine_weights_q2": {
+            "cosine": opt_weights_real_q2,
+            "sine": opt_weights_minus_imag_q2,
+        },
+        "opt_sine_weights_q2": {
+            "cosine": opt_weights_imag_q2,
+            "sine": opt_weights_real_q2,
+        },
+        "opt_minus_sine_weights_q2": {
+            "cosine": opt_weights_minus_imag_q2,
+            "sine": opt_weights_minus_real_q2,
+        },
+        "opt_cosine_weights_q1": {
+            "cosine": opt_weights_real_q1,
+            "sine": opt_weights_minus_imag_q1,
+        },
+        "opt_sine_weights_q1": {
+            "cosine": opt_weights_imag_q1,
+            "sine": opt_weights_real_q1,
+        },
+        "opt_minus_sine_weights_q1": {
+            "cosine": opt_weights_minus_imag_q1,
+            "sine": opt_weights_minus_real_q1,
+        },
+    },
+}
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/helper_functions.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/helper_functions.py
new file mode 100644
index 000000000..706f5a127
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/helper_functions.py	
@@ -0,0 +1,663 @@
+# helper functions pertinent to single- and two-qubit tomography
+
+import numpy as np
+import functools
+
+import matplotlib.pyplot as plt
+from matplotlib.colors import Normalize
+from matplotlib import cm, colormaps
+
+
+titlefont = {"color": "black", "weight": "normal", "size": 10}
+axisfont = {"color": "black", "weight": "normal", "size": 8}
+ticksfont = {"color": "black", "weight": "normal", "size": 8}
+
+
+from qm.qua import *
+
+
+def rotated_multiplexed_state_discrimination(I, I_st, Q, Q_st, states, states_st, resonators, thresholds):
+    """
+    Perform multiplexed state discrimination on two qubits
+
+    :param I: List of QUA variables to which the I quadrature measurements
+    are written
+    :param I_st: Optional list of QUA streams to which I values are streamed
+    :param Q: List of QUA variables to which the Q quadrature measurements
+    are written
+    :param states: List of QUA variables to which qubit state are written
+    :param states_st: Optional list of QUA streams to which qubit states
+    are streamed
+    :param resonators: List of integers denoting readout resonators
+    appearing in the config file as "ff{res}" elements
+    :param thresholds: List of readout thresholds from a calibrated
+    IQ blob two-state discriminator. I values above each threshold
+    signify that the qubit was found to be in the excited state
+
+    :return:
+    """
+
+    if type(resonators) is not list:
+        resonators = [resonators]
+
+    if type(thresholds) is not list:
+        thresholds = [thresholds]
+
+    for ind, res in enumerate(resonators):
+        measure(
+            "readout",
+            f"rr{res}",
+            None,
+            dual_demod.full("rotated_cos", "out1", "rotated_sin", "out2", I[ind]),
+            dual_demod.full("rotated_minus_sin", "out1", "rotated_cos", "out2", Q[ind]),
+        )
+
+        assign(states[ind], I[ind] > thresholds[ind])
+
+        if states_st is not None:
+            save(states[ind], states_st[ind])
+        if I_st is not None:
+            save(I[ind], I_st[ind])
+        if Q_st is not None:
+            save(Q[ind], Q_st[ind])
+
+
+@functools.lru_cache()
+def func_F1(n: int):
+    # one of the six gates to create one of the six cardinal Bloch sphere
+    # states
+    if n == 0:  # identity
+        return np.array([[1, 0], [0, 1]])
+    elif n == 1:  # X180/Y180
+        return np.array([[0, 1], [1, 0]])
+    elif n == 2:  # Y90
+        return (1 / np.sqrt(2)) * np.array([[1, -1], [1, 1]])
+    elif n == 3:  # Y-90
+        return (1 / np.sqrt(2)) * np.array([[1, 1], [-1, 1]])
+    elif n == 4:  # X-90
+        return (1 / np.sqrt(2)) * np.array([[1, 1j], [1j, 1]])
+    elif n == 5:  # X90
+        return (1 / np.sqrt(2)) * np.array([[1, -1j], [-1j, 1]])
+    else:  # Error
+        raise ValueError("Input integer should be between 0 and 5, inclusive")
+
+
+@functools.lru_cache()
+def func_E1(n: int):
+    # four Pauli operators
+    if n == 0:  # identity
+        return np.array([[1, 0], [0, 1]])
+    elif n == 1:  # X
+        return np.array([[0, 1], [1, 0]])
+    elif n == 2:  # Y
+        return np.array([[0, -1j], [1j, 0]])
+    elif n == 3:  # Z
+        return np.array([[1, 0], [0, -1]])
+    else:  # Error
+        raise ValueError("Input integer should be between 0 and 3, inclusive")
+
+
+@functools.lru_cache()
+def func_c1(i: int, j: int):
+    # Constants c1[i,j] such that
+    # func_E1[i] = sum_j c1[i,j]func_F1[j]|0><0|func_F1[j]^{\dagger}
+    if (i, j) == (0, 0):
+        return 1
+    elif (i, j) == (0, 1):
+        return 1
+    elif (i, j) == (1, 0):
+        return -(1 + 1j)
+    elif (i, j) == (1, 1):
+        return -(1 + 1j)
+    elif (i, j) == (1, 2):
+        return 2
+    elif (i, j) == (1, 4):
+        return 1j
+    elif (i, j) == (1, 5):
+        return 1j
+    elif (i, j) == (2, 4):
+        return 1
+    elif (i, j) == (2, 5):
+        return -1
+    elif (i, j) == (3, 0):
+        return 1
+    elif (i, j) == (3, 1):
+        return -1
+    else:
+        return 0
+
+
+def B_Bloch1(i, j, m, n):
+
+    # start from qubit ground state, create one of six Bloch sphere states
+    # using func_F1, add Pauli operator basis states func_E1, and measure one
+    # of the six Bloch sphere state projectors,
+    # func_F1^{\dagger} |0><0| func_F1. These constants appear when expressing
+    # the Pauli states in terms of combinations of Bloch sphere states,
+    # assuming you start from the qubit ground state
+
+    starting_state = np.array([[1, 0], [0, 0]])
+    bloch_state = func_F1(i) @ starting_state @ (func_F1(i).conj().T)
+    add_paulis = func_E1(m) @ bloch_state @ (func_E1(n).conj().T)
+    measure_bloch_state = (func_F1(j).conj().T) @ starting_state @ func_F1(j) @ add_paulis
+
+    return measure_bloch_state.trace()
+
+
+def P_Pauli1(l, k, m, n):
+
+    result = 0
+
+    for i in range(0, 6):
+        for j in range(0, 6):
+            result += func_c1(l, i) * np.conj(func_c1(k, j)) * B_Bloch1(i, j, m, n)
+
+    return result
+
+
+def map_from_bloch_state_to_pauli_basis1(l, k, arr):
+
+    if arr.shape != (6, 6):
+        raise ValueError("Input array must be 6x6")
+
+    if (k not in [0, 1, 2, 3]) or (l not in [0, 1, 2, 3]):
+        raise ValueError("Input indices must be between 0 and 3, inclusive")
+
+    result = 0
+
+    for i in range(0, 6):
+        for j in range(0, 6):
+            result += func_c1(l, i) * np.conj(func_c1(k, j)) * arr[i][j]
+
+    return result
+
+
+def plot_process_tomography1(chi_vector, save_file: str = None):
+
+    plt.rcParams["text.usetex"] = True
+
+    cmap = colormaps.get_cmap("viridis")
+
+    # set up figure
+    fig = plt.figure(figsize=(9, 3), dpi=250, facecolor="white")
+    ax = fig.add_subplot(111, projection="3d")
+
+    # coordinates
+    r = np.arange(4)
+    _x, _y = np.meshgrid(r, r)
+    x, y = _x.ravel(), _y.ravel()
+
+    values0 = np.abs(chi_vector)
+    values1 = np.angle(chi_vector)
+
+    top = values0
+    bottom = np.zeros_like(top)
+
+    width = depth = 0.7
+
+    norm = Normalize(vmin=-np.pi, vmax=np.pi)
+    colors = cmap(norm(values1))
+
+    xy_ticks_labels = [r"$I$", r"$X$", r"$Y$", r"$Z$"]
+
+    ax.bar3d(x, y, bottom, width, depth, top, shade=True, color=colors)
+    ax.set_xticks(r + 0.5, labels=xy_ticks_labels)
+    ax.set_yticks(r + 0.5, labels=xy_ticks_labels)
+    ax.set_zticks([0, (max(top) / 2).round(2), max(top).round(2)])
+    ax.set_xlabel("Prepared", fontdict=axisfont)
+    ax.set_ylabel("Measured", fontdict=axisfont)
+    ax.set_title(r"$\chi$ matrix", fontdict=titlefont)
+    ax.view_init(20, -60, 0)
+
+    sc = cm.ScalarMappable(cmap=cmap, norm=norm)
+    cbar = plt.colorbar(sc, ax=ax, pad=0.1, shrink=0.7)
+    cbar.set_ticks(
+        ticks=[-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi],
+        labels=[r"$-\pi$", r"$-\pi/2$", r"0", r"$\pi/2$", r"$\pi$"],
+    )
+
+    if save_file:
+        plt.savefig(save_file, bbox_inches="tight")
+
+    plt.show()
+
+
+@functools.lru_cache()
+def func_F2(n: int):
+    # one of the six gates to create one of the six cardinal Bloch sphere
+    # states
+    if n == 0:  # identity
+        return np.array([[1, 0], [0, 1]])
+    elif n == 1:  # X180/Y180
+        return np.array([[0, 1], [1, 0]])
+    elif n == 2:  # Y90
+        return (1 / np.sqrt(2)) * np.array([[1, -1], [1, 1]])
+    elif n == 3:  # Y-90
+        return (1 / np.sqrt(2)) * np.array([[1, 1], [-1, 1]])
+    elif n == 4:  # X-90
+        return (1 / np.sqrt(2)) * np.array([[1, 1j], [1j, 1]])
+    elif n == 5:  # X90
+        return (1 / np.sqrt(2)) * np.array([[1, -1j], [-1j, 1]])
+    else:  # Error
+        raise ValueError("Input integer should be between 0 and 5, inclusive")
+
+
+@functools.lru_cache()
+def func_E2(n: int):
+    # the 16 two-qubit Pauli operators
+
+    i = np.array([[1, 0], [0, 1]])
+    x = np.array([[0, 1], [1, 0]])
+    y = np.array([[0, -1j], [1j, 0]])
+    z = np.array([[1, 0], [0, -1]])
+
+    if n == 0:  # I x I
+        return np.kron(i, i)
+    elif n == 1:  # I x X
+        return np.kron(i, x)
+    elif n == 2:  # I x Y
+        return np.kron(i, y)
+    elif n == 3:  # I x Z
+        return np.kron(i, z)
+    elif n == 4:  # X x I
+        return np.kron(x, i)
+    elif n == 5:  # X x X
+        return np.kron(x, x)
+    elif n == 6:  # X x Y
+        return np.kron(x, y)
+    elif n == 7:  # X x Z
+        return np.kron(x, z)
+    elif n == 8:  # Y x I
+        return np.kron(y, i)
+    elif n == 9:  # Y x X
+        return np.kron(y, x)
+    elif n == 10:  # Y x Y
+        return np.kron(y, y)
+    elif n == 11:  # Y x Z
+        return np.kron(y, z)
+    elif n == 12:  # Z x I
+        return np.kron(z, i)
+    elif n == 13:  # Z x X
+        return np.kron(z, x)
+    elif n == 14:  # Z x Y
+        return np.kron(z, y)
+    elif n == 15:  # Z x Z
+        return np.kron(z, z)
+    else:
+        raise ValueError("Input integer should be between 0 and 15, inclusive")
+
+
+@functools.lru_cache()
+def func_c2(i: int, j: int, k: int):
+    # Constants c2[i, j, k] such that
+    # func_E2[i] = sum_{j,k} func_c2[i,j,k] (func_F2[j] x func_F2[k]) |0>|0><0|<0| (func_F2[j]^{\dagger} x func_F2[k]^{\dagger})
+    if (i, j, k) == (0, 0, 0):
+        return 1
+    elif (i, j, k) == (0, 0, 1):
+        return 1
+    elif (i, j, k) == (0, 1, 0):
+        return 1
+    elif (i, j, k) == (0, 1, 1):
+        return 1
+    elif (i, j, k) == (1, 0, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (1, 0, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (1, 0, 2):
+        return 2
+    elif (i, j, k) == (1, 0, 4):
+        return 1j
+    elif (i, j, k) == (1, 0, 5):
+        return 1j
+    elif (i, j, k) == (1, 1, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (1, 1, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (1, 1, 2):
+        return 2
+    elif (i, j, k) == (1, 1, 4):
+        return 1j
+    elif (i, j, k) == (1, 1, 5):
+        return 1j
+    elif (i, j, k) == (2, 0, 4):
+        return 1
+    elif (i, j, k) == (2, 0, 5):
+        return -1
+    elif (i, j, k) == (2, 1, 4):
+        return 1
+    elif (i, j, k) == (2, 1, 5):
+        return -1
+    elif (i, j, k) == (3, 0, 0):
+        return 1
+    elif (i, j, k) == (3, 0, 1):
+        return -1
+    elif (i, j, k) == (3, 1, 0):
+        return 1
+    elif (i, j, k) == (3, 1, 1):
+        return -1
+    elif (i, j, k) == (4, 0, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (4, 0, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (4, 1, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (4, 1, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (4, 2, 0):
+        return 2
+    elif (i, j, k) == (4, 2, 1):
+        return 2
+    elif (i, j, k) == (4, 4, 0):
+        return 1j
+    elif (i, j, k) == (4, 4, 1):
+        return 1j
+    elif (i, j, k) == (4, 5, 0):
+        return 1j
+    elif (i, j, k) == (4, 5, 1):
+        return 1j
+    elif (i, j, k) == (5, 0, 0):
+        return (1 + 1j) ** 2
+    elif (i, j, k) == (5, 0, 1):
+        return (1 + 1j) ** 2
+    elif (i, j, k) == (5, 1, 0):
+        return (1 + 1j) ** 2
+    elif (i, j, k) == (5, 1, 1):
+        return (1 + 1j) ** 2
+    elif (i, j, k) == (5, 0, 2):
+        return -2 * (1 + 1j)
+    elif (i, j, k) == (5, 1, 2):
+        return -2 * (1 + 1j)
+    elif (i, j, k) == (5, 2, 0):
+        return -2 * (1 + 1j)
+    elif (i, j, k) == (5, 2, 1):
+        return -2 * (1 + 1j)
+    elif (i, j, k) == (5, 0, 4):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 0, 5):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 1, 4):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 1, 5):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 4, 0):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 4, 1):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 5, 0):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 5, 1):
+        return -1j * (1 + 1j)
+    elif (i, j, k) == (5, 2, 2):
+        return 4
+    elif (i, j, k) == (5, 2, 4):
+        return 2 * 1j
+    elif (i, j, k) == (5, 2, 5):
+        return 2 * 1j
+    elif (i, j, k) == (5, 4, 2):
+        return 2 * 1j
+    elif (i, j, k) == (5, 5, 2):
+        return 2 * 1j
+    elif (i, j, k) == (5, 4, 4):
+        return -1
+    elif (i, j, k) == (5, 4, 5):
+        return -1
+    elif (i, j, k) == (5, 5, 4):
+        return -1
+    elif (i, j, k) == (5, 5, 5):
+        return -1
+    elif (i, j, k) == (6, 0, 4):
+        return -(1 + 1j)
+    elif (i, j, k) == (6, 1, 4):
+        return -(1 + 1j)
+    elif (i, j, k) == (6, 0, 5):
+        return 1 + 1j
+    elif (i, j, k) == (6, 1, 5):
+        return 1 + 1j
+    elif (i, j, k) == (6, 2, 4):
+        return 2
+    elif (i, j, k) == (6, 2, 5):
+        return -2
+    elif (i, j, k) == (6, 4, 4):
+        return 1j
+    elif (i, j, k) == (6, 5, 4):
+        return 1j
+    elif (i, j, k) == (6, 4, 5):
+        return -1j
+    elif (i, j, k) == (6, 5, 5):
+        return -1j
+    elif (i, j, k) == (7, 0, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (7, 1, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (7, 0, 1):
+        return 1 + 1j
+    elif (i, j, k) == (7, 1, 1):
+        return 1 + 1j
+    elif (i, j, k) == (7, 2, 0):
+        return 2
+    elif (i, j, k) == (7, 2, 1):
+        return -2
+    elif (i, j, k) == (7, 4, 0):
+        return 1j
+    elif (i, j, k) == (7, 5, 0):
+        return 1j
+    elif (i, j, k) == (7, 4, 1):
+        return -1j
+    elif (i, j, k) == (7, 5, 1):
+        return -1j
+    elif (i, j, k) == (8, 4, 0):
+        return 1
+    elif (i, j, k) == (8, 4, 1):
+        return 1
+    elif (i, j, k) == (8, 5, 0):
+        return -1
+    elif (i, j, k) == (8, 5, 1):
+        return -1
+    elif (i, j, k) == (9, 4, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (9, 4, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (9, 5, 0):
+        return 1 + 1j
+    elif (i, j, k) == (9, 5, 1):
+        return 1 + 1j
+    elif (i, j, k) == (9, 4, 4):
+        return 1j
+    elif (i, j, k) == (9, 4, 5):
+        return 1j
+    elif (i, j, k) == (9, 5, 4):
+        return -1j
+    elif (i, j, k) == (9, 5, 5):
+        return -1j
+    elif (i, j, k) == (9, 4, 2):
+        return 2
+    elif (i, j, k) == (9, 5, 2):
+        return -2
+    elif (i, j, k) == (10, 4, 4):
+        return 1
+    elif (i, j, k) == (10, 5, 5):
+        return 1
+    elif (i, j, k) == (10, 4, 5):
+        return -1
+    elif (i, j, k) == (10, 5, 4):
+        return -1
+    elif (i, j, k) == (11, 4, 0):
+        return 1
+    elif (i, j, k) == (11, 5, 1):
+        return 1
+    elif (i, j, k) == (11, 4, 1):
+        return -1
+    elif (i, j, k) == (11, 5, 0):
+        return -1
+    elif (i, j, k) == (12, 0, 0):
+        return 1
+    elif (i, j, k) == (12, 0, 1):
+        return 1
+    elif (i, j, k) == (12, 1, 0):
+        return -1
+    elif (i, j, k) == (12, 1, 1):
+        return -1
+    elif (i, j, k) == (13, 0, 0):
+        return -(1 + 1j)
+    elif (i, j, k) == (13, 0, 1):
+        return -(1 + 1j)
+    elif (i, j, k) == (13, 1, 0):
+        return 1 + 1j
+    elif (i, j, k) == (13, 1, 1):
+        return 1 + 1j
+    elif (i, j, k) == (13, 0, 2):
+        return 2
+    elif (i, j, k) == (13, 1, 2):
+        return -2
+    elif (i, j, k) == (13, 0, 4):
+        return 1j
+    elif (i, j, k) == (13, 0, 5):
+        return 1j
+    elif (i, j, k) == (13, 1, 4):
+        return -1j
+    elif (i, j, k) == (13, 1, 5):
+        return -1j
+    elif (i, j, k) == (14, 0, 4):
+        return 1
+    elif (i, j, k) == (14, 1, 5):
+        return 1
+    elif (i, j, k) == (14, 0, 5):
+        return -1
+    elif (i, j, k) == (14, 1, 4):
+        return -1
+    elif (i, j, k) == (15, 0, 0):
+        return 1
+    elif (i, j, k) == (15, 1, 1):
+        return 1
+    elif (i, j, k) == (15, 0, 1):
+        return -1
+    elif (i, j, k) == (15, 1, 0):
+        return -1
+    else:
+        return 0
+
+
+def B_Bloch2(i, j, k, l, m, n):
+
+    # start from qubits ground states, create one of six Bloch sphere states
+    # using func_F2 for each of the two qubits, add one of the 16
+    # Pauli operator basis func_E2, and measure one of the six Bloch sphere
+    # state projectors for each qubit, i.e.
+    # func_F2[.]^{\dagger}func_F2[.]^{\dagger} |0>|0><0|<0| func_F2[.] func_F2[.].
+    # These B_Bloch2 constants appear when expressing the Pauli states
+    # in terms of combinations of Bloch sphere states, assuming you
+    # start from the qubits' ground states
+
+    starting_state = np.array([[1, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])
+
+    bloch_state = np.kron(func_F2(i), func_F2(j)) @ starting_state @ np.kron(func_F2(i).conj().T, func_F2(j).conj().T)
+    add_paulis = func_E2(m) @ bloch_state @ (func_E2(n).conj().T)
+    measure_bloch_state = (
+        np.kron(func_F2(k).conj().T, func_F2(l).conj().T)
+        @ starting_state
+        @ np.kron(func_F2(k), func_F2(l))
+        @ add_paulis
+    )
+
+    return measure_bloch_state.trace()
+
+
+def P_Pauli2(s, t, m, n):
+
+    result = 0
+
+    for i in range(0, 6):
+        for j in range(0, 6):
+            for k in range(0, 6):
+                for l in range(0, 6):
+                    result += func_c2(s, i, j) * np.conj(func_c2(t, k, l)) * B_Bloch2(i, j, k, l, m, n)
+
+    return result
+
+
+def map_from_bloch_state_to_pauli2(q, n, arr):
+
+    if arr.shape != (6, 6, 6, 6):
+        raise ValueError("Input array must be 6x6x6x6")
+
+    if (q not in range(0, 16)) or (n not in range(0, 16)):
+        raise ValueError("Input indices must be between 0 and 15, inclusive")
+
+    result = 0
+
+    for i in range(0, 6):  # first qubit prepare
+        for j in range(0, 6):  # second qubit prepare
+            for k in range(0, 6):  # first qubit measure
+                for l in range(0, 6):  # second qubit measure
+                    result += func_c2(q, i, j) * np.conj(func_c2(n, k, l)) * arr[i][j][k][l]
+
+    return result
+
+
+def plot_process_tomography2(chi_vector, save_file: str = None):
+
+    plt.rcParams["text.usetex"] = True
+
+    ticksfont = {"color": "black", "weight": "normal", "size": 4}
+
+    cmap = colormaps.get_cmap("viridis")
+
+    # set up figure
+    fig = plt.figure(figsize=(9, 3), dpi=250, facecolor="white")
+    ax = fig.add_subplot(111, projection="3d")
+
+    # coordinates
+    r = np.arange(16)
+    _x, _y = np.meshgrid(r, r)
+    x, y = _x.ravel(), _y.ravel()
+
+    values0 = np.abs(chi_vector)
+    values1 = np.angle(chi_vector)
+
+    top = values0
+    bottom = np.zeros_like(top)
+
+    width = depth = 0.7
+
+    norm = Normalize(vmin=-np.pi, vmax=np.pi)
+    colors = cmap(norm(values1))
+
+    xy_ticks_labels = [
+        r"$II$",
+        r"$IX$",
+        r"$IY$",
+        r"$IZ$",
+        r"$XI$",
+        r"$XX$",
+        r"$XY$",
+        r"$XZ$",
+        r"$YI$",
+        r"$YX$",
+        r"$YY$",
+        r"$YZ$",
+        r"$ZI$",
+        r"$ZX$",
+        r"$ZY$",
+        r"$ZZ$",
+    ]
+
+    ax.bar3d(x, y, bottom, width, depth, top, shade=True, color=colors)
+    ax.set_xticks(r + 0.5, labels=xy_ticks_labels, fontdict=ticksfont)
+    ax.set_yticks(r + 0.5, labels=xy_ticks_labels, fontdict=ticksfont)
+    ax.set_zticks([0, (max(top) / 2).round(2), max(top).round(2)])
+    ax.set_xlabel("Prepared", fontdict=axisfont)
+    ax.set_ylabel("Measured", fontdict=axisfont)
+    ax.set_title(r"$\chi$ matrix", fontdict=titlefont)
+    ax.view_init(20, -60, 0)
+
+    sc = cm.ScalarMappable(cmap=cmap, norm=norm)
+    cbar = plt.colorbar(sc, ax=ax, pad=0.1, shrink=0.7)
+    cbar.set_ticks(
+        ticks=[-np.pi, -np.pi / 2, 0, np.pi / 2, np.pi],
+        labels=[r"$-\pi$", r"$-\pi/2$", r"0", r"$\pi/2$", r"$\pi$"],
+    )
+
+    if save_file:
+        plt.savefig(save_file, bbox_inches="tight")
+
+    plt.show()
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_ideal.png b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_ideal.png
new file mode 100644
index 000000000..40bac7ccb
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_ideal.png differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_real.png b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_real.png
new file mode 100644
index 000000000..e8f217ba0
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/iswap_real.png differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/setup.png b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/setup.png
new file mode 100644
index 000000000..86f7b94ac
Binary files /dev/null and b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/setup.png differ
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-process-tomography.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-process-tomography.py
new file mode 100644
index 000000000..629ca0e96
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-process-tomography.py	
@@ -0,0 +1,242 @@
+#!/usr/bin/env python
+
+"""
+        SINGLE QUBIT PROCESS TOMOGRAPHY
+The sequence consists of preparing the qubit into one of the six cardinal Bloch sphere
+states using calibrated gates from the config, applying the operation/process - typically
+ideally a unitary matrix - under investigation, and then measuring the state of the qubit,
+by way of the readout resonator, in the X, Y and Z bases (specifically, measuring the
+projector of one of the same six cardinal Bloch sphere states). This is repeated for
+the full set of input states and measurement projectors, which is tomographically
+complete. The output, 'probs', containing the 'Bloch sphere' prepared states and
+measurement projectors, is mapped to the Pauli basis (giving a 'measurement vector'),
+and the equation "measurement_vector = B \times chi_vector" (where B is a constant)
+inverted to give the chi process matrix for the process under
+investigation.
+
+
+
+Prerequisites:
+    - Having found the resonance frequency of the resonator coupled to the qubit under study (resonator_spectroscopy)
+    - Having calibrated qubit pi and pi/2 pulses by running qubit spectroscopy, rabi_chevron, power_rabi and updating the config
+    - Having calibrated the readout (readout_frequency, amplitude, duration_optimization IQ_blobs) for better SNR, and having
+      saved the derived readout threshold value in the config
+    - Set the desired flux bias in the case of flux-tunable qubits
+
+This script implements the logic outlined in the notebook
+https://github.com/bornman-nick/quantum-state-and-process-tomography. See
+that notebook for a detailed explanation of standard single qubit process
+tomography.
+
+Note that this routine carries out standard process tomography on a single qubit.
+Given that the measurements recorded are simply whether the final state of the
+qubit is in the ground state or not for each state preparation and measurement
+setting (instead of also deciding whether the qubit is perhaps in the excited
+state and incrementing the counter of a complementary measurement setting),
+this routine scales as 6**2, rather than 6*3, as is the case in
+most formulations of standard process tomography. However, the routine below
+is fast enough that this doesn't matter too much.
+"""
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool
+import numpy as np
+from scipy.linalg import solve
+
+from helper_functions import (
+    P_Pauli1,
+    plot_process_tomography1,
+    map_from_bloch_state_to_pauli_basis1,
+    func_E1,
+)
+
+
+###################
+# The QUA program #
+###################
+
+# qubit under test, assuming there are multiple qubits
+# with XY line elements "q<qubit>_xy", flux lines
+# "q<qubit>_z", and readout resonator elements "rr<qubit>"
+qubit = 1
+
+if qubit == 1:
+    threshold = ge_threshold_q1
+elif qubit == 2:
+    threshold = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit number chosen")
+
+
+n_avg = 10_000
+
+
+# subroutine to prepare desired qubit gate/process
+def analysed_process(qubit):
+    # write whatever QUA code you need, here, in order to perform the
+    # desired process which we want to subject to tomography.
+    # For example, to analyse the Y gate:
+    play("y180", f"q{qubit}_xy")
+
+
+ideal_gate = func_E1(2)
+
+
+# subroutine to switch between preparing the tomographically-complete
+# set of input states
+def prepare_state(i):
+    with switch_(i):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit}_xy")
+        with case_(1):
+            play("x180", f"q{qubit}_xy")
+        with case_(2):
+            play("y90", f"q{qubit}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit}_xy")
+        with case_(5):
+            play("x90", f"q{qubit}_xy")
+
+
+# subroutine to switch between the tomographically-complete
+# set of bases in which to measure
+def measurement_basis_change(j):
+    with switch_(j):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit}_xy")
+        with case_(1):
+            play("x180", f"q{qubit}_xy")
+        with case_(2):
+            play("y90", f"q{qubit}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit}_xy")
+        with case_(5):
+            play("x90", f"q{qubit}_xy")
+
+
+with program() as single_qubit_process_tomography:
+
+    n = declare(int)
+    n_st = declare_stream()
+
+    I = declare(fixed)
+    Q = declare(fixed)
+
+    state = declare(bool)  # QUA variable for the measured qubit state
+    state_st = declare_stream()  # Stream for the qubit state
+
+    c = declare(int)  # QUA variable for switching between state preparation/creations
+    m = declare(int)  # QUA variable for switching between basis projections/measurements
+
+    with for_(n, 0, n < n_avg, n + 1):  # QUA for_ loop for averaging
+        with for_(c, 0, c <= 5, c + 1):  # QUA for_ loop for switching between state preparations
+            with for_(m, 0, m <= 5, m + 1):  # QUA for_ loop for switching between basis projections/measurements
+
+                # prepare qubit in one of six cardinal Bloch sphere states
+                prepare_state(c)
+
+                align()
+
+                # the process to be analysed
+                analysed_process(qubit)
+
+                align()
+
+                # projective measurement basis change
+                measurement_basis_change(m)
+
+                align()
+
+                measure(
+                    "readout",
+                    f"rr{qubit}",
+                    None,
+                    dual_demod.full("rotated_cos", "out1", "rotated_sin", "out2", I),
+                    dual_demod.full("rotated_minus_sin", "out1", "rotated_cos", "out2", Q),
+                )
+
+                align()
+
+                # track number of times the qubit is found to be
+                # in the ground state
+                assign(state, I < threshold)
+
+                save(state, state_st)
+
+                wait(thermalization_time * u.ns)
+
+        save(n, n_st)
+
+    with stream_processing():
+        n_st.save("iteration")
+        state_st.boolean_to_int().buffer(6).buffer(6).average().save("probs")
+
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+
+simulate = False
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=10_000)  # In clock cycles = 4ns
+    job = qmm.simulate(config, single_qubit_process_tomography, simulation_config)
+    job.get_simulated_samples().con1.plot()
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(single_qubit_process_tomography)
+    # Get results from QUA program
+    results = fetching_tool(job, data_list=["iteration", "probs"], mode="live")
+
+    while results.is_processing():
+
+        # Fetch results
+        iteration, probs = results.fetch_all()
+
+        # Progress bar
+        progress_counter(iteration, n_avg, start_time=results.get_start_time())
+
+    # Close the quantum machines at the end in order to put all flux biases to 0 so that the fridge doesn't heat-up
+    qm.close()
+
+    # post-processing
+    PMatrix = np.array(
+        [
+            [
+                P_Pauli1(
+                    np.floor(v / 4).astype(int),
+                    v % 4,
+                    np.floor(w / 4).astype(int),
+                    w % 4,
+                )
+                for w in range(16)
+            ]
+            for v in range(16)
+        ]
+    )
+
+    pauli_basis_measurements = lambda l, k: map_from_bloch_state_to_pauli_basis1(l, k, probs)
+
+    measurement_vector = np.array([pauli_basis_measurements(np.floor(v / 4).astype(int), v % 4) for v in range(16)])
+
+    chi_vector = solve(PMatrix, measurement_vector)
+    chi_matrix = chi_vector.reshape(4, 4)
+
+    plot_process_tomography1(chi_vector)
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-state-tomography.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-state-tomography.py
new file mode 100644
index 000000000..dd9ff87e8
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/single-qubit-state-tomography.py	
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+
+"""
+        SINGLE QUBIT STATE TOMOGRAPHY
+The sequence consists of preparing the qubit into a chosen state using calibrated gates from the config, and
+measuring the state of the qubit, by way of the readout resonator, in the X, Y and Z bases. The output,
+'probs', is a triple of probability 'difference' counts for measuring the qubit in one eigenstate of the
+X/Y/Z bases, minus the probability of measuring it in that eigenstate's complement.
+
+These values are scaled into the usual Stokes parameters and used to infer the state of the qubit; see
+https://research.physics.illinois.edu/QI/Photonics/tomography-files/tomo_chapter_2004.pdf for further details
+
+Note that this program is similar to 
+qua-libs/Quantum-Control-Applications/Superconducting/Single-Fixed-Transmon
+/19_state_tomography.py, which the author became aware of after having
+written the current script
+
+Prerequisites:
+    - Having found the resonance frequency of the resonator coupled to the qubit under study (resonator_spectroscopy)
+    - Having calibrated qubit pi and pi/2 pulses by running qubit spectroscopy, rabi_chevron, power_rabi and updating the config
+    - Having calibrated the readout (readout_frequency, amplitude, duration_optimization IQ_blobs) for better SNR, and having
+      saved the derived readout threshold value in the config
+    - Set the desired flux bias in the case of flux-tunable qubits
+
+This script implements the logic outlined in the notebook
+https://github.com/bornman-nick/quantum-state-and-process-tomography. See that
+notebook for a detailed explanation of standard single qubit state
+tomography.
+"""
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool
+import numpy as np
+
+
+###################
+# The QUA program #
+###################
+
+# qubit under test, assuming there are multiple qubits
+# with XY line elements "q<qubit>_xy", flux lines
+# "q<qubit>_z", and readout resonator elements "rr<qubit>"
+qubit = 1
+
+if qubit == 1:
+    threshold = ge_threshold_q1
+elif qubit == 2:
+    threshold = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit number chosen")
+
+
+n_avg = 10_000
+
+
+# subroutine to prepare desired qubit state
+def prepare_state(qubit):
+    # write whatever QUA code you need in order to create the
+    # state to perform tomography on, from an initial
+    # ground state. For example, to create the |1> state
+    play("y180", f"q{qubit}_xy")
+
+
+with program() as single_qubit_state_tomography:
+    n = declare(int)  # QUA variable for average loop
+    n_st = declare_stream()  # Stream for the averaging iteration 'n'
+
+    state = declare(bool)  # QUA variable for the qubit state
+    state_st = declare_stream()  # Stream for the qubit state
+
+    p = declare(int)  # QUA variable for switching between projections
+
+    I = declare(fixed)
+    Q = declare(fixed)
+
+    with for_(n, 0, n < n_avg, n + 1):  # QUA for_ loop for averaging
+        with for_(p, 0, p <= 2, p + 1):  # QUA for_ loop for switching between basis changes
+
+            prepare_state(qubit)
+            align()
+
+            with switch_(c):
+                with case_(0):  # basis X
+
+                    # Map the X-component of the Bloch vector onto the Z-axis
+                    # 1/sqrt(2)(|0>+|1>) -> |0>; 1/sqrt(2)(|0>-|1>) -> |1>
+                    play("-y90", f"q{qubit}_xy")
+
+                    align(f"q{qubit}_xy", f"rr{qubit}")
+
+                with case_(1):  # basis Y
+
+                    # Map the Y-component of the Bloch vector onto the Z-axis
+                    # 1/sqrt(2)(|0>+i|1>) -> |0>; 1/sqrt(2)(|0>-i|1>) -> |1>
+                    play("x90", f"q{qubit}_xy")
+
+                    align(f"q{qubit}_xy", f"rr{qubit}")
+
+                with case_(2):  # basis Z
+
+                    align(f"q{qubit}_xy", f"rr{qubit}")
+
+                measure(
+                    "readout",
+                    f"rr{qubit}",
+                    None,
+                    dual_demod.full("rotated_cos", "out1", "rotated_sin", "out2", I),
+                    dual_demod.full("rotated_minus_sin", "out1", "rotated_cos", "out2", Q),
+                )
+
+            # True if qubit state is |1>, False if |0>
+
+            assign(state, I > threshold)
+
+            wait(thermalization_time * u.ns)
+
+            save(state, state_st)
+
+        save(n, n_st)
+
+    with stream_processing():
+        n_st.save("iteration")
+        state_st.boolean_to_int().buffer(3).average().save("probs")
+
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+simulate = False
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=10_000)  # In clock cycles = 4ns
+    job = qmm.simulate(config, single_qubit_state_tomography, simulation_config)
+    job.get_simulated_samples().con1.plot()
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(single_qubit_state_tomography)
+    # Get results from QUA program
+    results = fetching_tool(job, data_list=["probs", "iteration"], mode="live")
+
+    while results.is_processing():
+        # Fetch results
+        probs, iteration = results.fetch_all()
+        # Progress bar
+        progress_counter(iteration, n_avg, start_time=results.get_start_time())
+
+        # Converts the (0,1) -> |g>,|e> convention, arising from the I>threshold
+        # assignment, to (1,-1) -> |g>,|e>, which aligns with the Stokes parameter
+        # definitions from the projector probabilities
+        prob = -2 * (probs - 0.5)
+
+    # Close the quantum machines at the end in order to put all flux biases to 0 so that the fridge doesn't heat-up
+    qm.close()
+
+    # Reconstruct the density matrix
+    I = np.array([[1, 0], [0, 1]])
+    sigma_x = np.array([[0, 1], [1, 0]])
+    sigma_y = np.array([[0, -1j], [1j, 0]])
+    sigma_z = np.array([[1, 0], [0, -1]])
+
+    rho = 0.5 * (I + prob[0] * sigma_x + prob[1] * sigma_y + prob[2] * sigma_z)
+    print(f"The density matrix is:\n{rho}")
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-process-tomography.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-process-tomography.py
new file mode 100644
index 000000000..b83509a70
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-process-tomography.py	
@@ -0,0 +1,331 @@
+#!/usr/bin/env python
+
+"""
+        TWO QUBIT PROCESS TOMOGRAPHY
+The sequence consists of preparing the qubits each into one of their six cardinal Bloch sphere
+states using calibrated gates from the config, applying the operation/process - typically
+ideally a unitary matrix - under investigation, and then measuring the state of the qubits,
+by way of the readout resonator, in their X, Y and Z bases (specifically, measuring the
+projector of one of the same six cardinal Bloch sphere states). This is repeated for
+the full set of input states and measurement projectors, which is tomographically
+complete. The output, 'probs', containing the 'Bloch sphere' prepared states and
+measurement projectors, is mapped to the Pauli basis (giving a 'measurement vector'),
+and the equation "measurement_vector = B \times chi_vector" (where B is a constant)
+inverted to give the chi process matrix for the process under
+investigation.
+
+Computing PMatrix from scratch in the two-qubit process case takes a while,
+so included with this project is a serialised file of Pmatrix
+
+Prerequisites:
+    - Having found the resonance frequency of the resonator coupled to the qubits under study (resonator_spectroscopy)
+    - Having calibrated qubits' pi and pi/2 pulses by running qubit spectroscopy, rabi_chevron, power_rabi and updating the config
+    - Having calibrated the readout (readout_frequency, amplitude, duration_optimization IQ_blobs) for better SNR, and having
+      saved the derived readout threshold values in the config
+    - Set the desired flux biases in the case of flux-tunable qubits
+
+This script implements the logic outlined in the notebook
+https://github.com/bornman-nick/quantum-state-and-process-tomography. See that
+notebook for a detailed explanation of standard two qubit process tomography
+
+Note that this routine carries out standard process tomography for two qubits.
+Given that the measurements recorded are simply whether the final state of the
+qubits are in their ground state or not for each state preparation and measurement
+setting (instead of also deciding whether the qubits are perhaps in the excited
+states and incrementing the counter of a complementary measurement setting),
+this routine scales as 6**4, rather than (6**2)*(3**2), as is the case in
+most formulations of standard process tomography. However, the routine below
+is fast enough that this doesn't matter too much.
+"""
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool
+import numpy as np
+from scipy.linalg import solve
+
+import os
+import pickle
+
+from helper_functions import (
+    P_Pauli2,
+    plot_process_tomography2,
+    map_from_bloch_state_to_pauli_basis2,
+    func_F2,
+)
+
+
+###################
+# The QUA program #
+###################
+
+# load PMatrix if pickle file is present
+# note: this takes a while to compute if you do not load the serialised
+# PMatrix file
+file_location = "<path to PMatrix2 file>"  # ./PMatrix2.pkl
+
+if os.path.isfile(file_location):
+    with open(file_location, "rb") as file:
+        PMatrix = pickle.load(file)
+else:
+    # computing this takes a while
+    PMatrix = np.array(
+        [
+            [
+                P_Pauli2(
+                    np.floor(v / 16).astype(int),
+                    v % 16,
+                    np.floor(w / 16).astype(int),
+                    w % 16,
+                )
+                for w in range(16**2)
+            ]
+            for v in range(16**2)
+        ]
+    )
+
+# qubits under test, assuming there are multiple qubits
+# with XY line elements "q<qubit>_xy", flux lines
+# "q<qubit>_z", and readout resonator elements "rr<qubit>"
+qubit1 = 1
+qubit2 = 2
+
+if qubit1 == 1:
+    threshold1 = ge_threshold_q1
+elif qubit1 == 2:
+    threshold1 = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit1 number chosen")
+
+if qubit2 == 1:
+    threshold2 = ge_threshold_q1
+elif qubit2 == 2:
+    threshold2 = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit2 number chosen")
+
+if qubit1 == qubit2:
+    raise ValueError(f"The value of qubit1 cannot equal that of qubit2")
+
+
+n_avg = 5_000
+
+
+# subroutine to prepare desired gate/process
+def analysed_process(qubit1, qubit2):
+    # write whatever QUA code you need, here, in order to perform the
+    # desired process which we want to subject to tomography.
+    # For example, to analyse the X gate on qubit1 and the -Y90
+    # gate on qubit2:
+    play("x180", f"q{qubit1}_xy")
+    play("-y90", f"q{qubit2}_xy")
+
+
+ideal_gate = np.kron(func_F2(1), func_F2(3))
+
+
+# subroutine to switch between preparing the tomographically-complete
+# set of input states
+def prepare_states(i, j, qubit1, qubit2):
+
+    # Prepare Bloch state i on qubit1
+    with switch_(i):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit1}_xy")
+        with case_(1):
+            play("x180", f"q{qubit1}_xy")
+        with case_(2):
+            play("y90", f"q{qubit1}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit1}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit1}_xy")
+        with case_(5):
+            play("x90", f"q{qubit1}_xy")
+
+    # Prepare Bloch state j on qubit2
+    with switch_(j):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit2}_xy")
+        with case_(1):
+            play("x180", f"q{qubit2}_xy")
+        with case_(2):
+            play("y90", f"q{qubit2}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit2}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit2}_xy")
+        with case_(5):
+            play("x90", f"q{qubit2}_xy")
+
+
+# subroutine to switch between the tomographically-complete
+# set of bases in which to measure
+def measurement_basis_change(k, l, qubit1, qubit2):
+
+    # Change basis with operation k on qubit1
+    with switch_(k):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit2}_xy")
+        with case_(1):
+            play("x180", f"q{qubit1}_xy")
+        with case_(2):
+            play("y90", f"q{qubit1}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit1}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit1}_xy")
+        with case_(5):
+            play("x90", f"q{qubit1}_xy")
+
+    # Change basis with operation l on qubit2
+    with switch_(l):
+        with case_(0):
+            wait(pi_len // 4, f"q{qubit2}_xy")
+        with case_(1):
+            play("x180", f"q{qubit2}_xy")
+        with case_(2):
+            play("y90", f"q{qubit2}_xy")
+        with case_(3):
+            play("-y90", f"q{qubit2}_xy")
+        with case_(4):
+            play("-x90", f"q{qubit2}_xy")
+        with case_(5):
+            play("x90", f"q{qubit2}_xy")
+
+
+with program() as two_qubit_process_tomography:
+
+    n = declare(int)
+    n_st = declare_stream()
+
+    I1 = declare(fixed)
+    Q1 = declare(fixed)
+    I2 = declare(fixed)
+    Q2 = declare(fixed)
+
+    state = declare(bool)  # QUA variable for the measured qubits state
+    state_st = declare_stream()  # Stream for the qubits state
+
+    c1 = declare(int)  # QUA variable for switching between state preparation/creations on qubit 1
+    c2 = declare(int)  # QUA variable for switching between state preparation/creations on qubit 2
+    m1 = declare(int)  # QUA variable for switching between Bloch basis projections/measurements on qubit 1
+    m2 = declare(int)  # QUA variable for switching between Bloch basis projections/measurements on qubit 2
+
+    with for_(n, 0, n < n_avg, n + 1):  # QUA for_ loop for averaging
+        with for_(c1, 0, c1 <= 5, c1 + 1):  # QUA for_ loop for switching between state preparations on qubit 1
+            with for_(c2, 0, c2 <= 5, c2 + 1):  # QUA for_ loop for switching between state preparations on qubit 2
+                with for_(
+                    m1, 0, m1 <= 5, m1 + 1
+                ):  # QUA for_ loop for switching between Bloch basis projections/measurements on qubit 1
+                    with for_(
+                        m2, 0, m2 <= 5, m2 + 1
+                    ):  # QUA for_ loop for switching between Bloch basis projections/measurements on qubit 2
+
+                        reset_frame(f"q{qubit1}_xy")
+                        reset_frame(f"q{qubit2}_xy")
+
+                        reset_phase(f"q{qubit1}_xy")
+                        reset_phase(f"q{qubit2}_xy")
+
+                        # prepare qubit1 and qubit 2 in one of six Bloch sphere states
+                        prepare_states(c1, c2, qubit1, qubit2)
+
+                        align()
+
+                        # apply the process to be analysed
+                        analysed_process(qubit1, qubit2)
+
+                        aling()
+
+                        # projective measurement basis change
+                        measurement_basis_change(m1, m2, qubit1, qubit2)
+
+                        align()
+
+                        measure(
+                            "readout",
+                            f"rr{qubit1}",
+                            None,
+                            dual_demod.full("rotated_cos", "out1", "rotated_sin", "out2", I1),
+                            dual_demod.full("rotated_minus_sin", "out1", "rotated_cos", "out2", Q1),
+                        )
+
+                        measure(
+                            "readout",
+                            f"rr{qubit2}",
+                            None,
+                            dual_demod.full("rotated_cos", "out1", "rotated_sin", "out2", I2),
+                            dual_demod.full("rotated_minus_sin", "out1", "rotated_cos", "out2", Q2),
+                        )
+
+                        align()
+
+                        # track the number of clicks when both qubit 1 and qubit 2 are
+                        # in their ground states
+                        with if_((I1 < threshold1) & (I2 < threshold2)):
+                            assign(state, True)
+                        with else_():
+                            assign(state, False)
+
+                        save(state, state_st)
+
+                        wait(thermalization_time * u.ns, f"rr{qubit1}", f"rr{qubit2}")
+
+        save(n, n_st)
+
+    with stream_processing():
+        n_st.save("iteration")
+        state_st.boolean_to_int().buffer(6).buffer(6).buffer(6).buffer(6).average().save("probs")
+
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+
+simulate = False
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=10_000)  # In clock cycles = 4ns
+    job = qmm.simulate(config, two_qubit_process_tomography, simulation_config)
+    job.get_simulated_samples().con1.plot()
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(two_qubit_process_tomography)
+    # Get results from QUA program
+    results = fetching_tool(job, data_list=["iteration", "probs"], mode="live")
+
+    while results.is_processing():
+
+        # Fetch results
+        iteration, probs = results.fetch_all()
+
+        # Progress bar
+        progress_counter(iteration, n_avg, start_time=results.get_start_time())
+
+    # Close the quantum machines at the end in order to put all flux biases to 0 so that the fridge doesn't heat-up
+    qm.close()
+
+    # post-processing
+    pauli_basis_measurements = lambda q, n: map_from_bloch_state_to_pauli_basis2(q, n, probs)
+
+    measurement_vector = np.array(
+        [pauli_basis_measurements(np.floor(v / 16).astype(int), v % 16) for v in range(16**2)]
+    )
+
+    chi_vector = solve(PMatrix, measurement_vector)
+    chi_matrix = chi_vector.reshape(16, 16)
+
+    plot_process_tomography2(chi_vector)
diff --git a/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-state-tomography.py b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-state-tomography.py
new file mode 100644
index 000000000..f6d0c17eb
--- /dev/null
+++ b/Quantum-Control-Applications/Superconducting/Two-Flux-Tunable-Coupled-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/two-qubit-state-tomography.py	
@@ -0,0 +1,349 @@
+#!/usr/bin/env python
+
+"""
+        TWO QUBIT STATE TOMOGRAPHY
+The sequence consists of preparing the qubits into a chosen state using
+calibrated gates from the config, and measuring the state of the qubits, by way
+of the readout resonator, in their X, Y and Z bases. The output, 'probs', is a
+list of 15 probability 'difference' counts corresponding with the 15 Stokes
+parameters for two qubit state tomography. These values are used to infer the
+two-qubit state; see
+https://research.physics.illinois.edu/QI/Photonics/tomography-files/tomo_chapter_2004.pdf
+for further details.
+
+Prerequisites:
+    - Having found the resonance frequencies of the resonators coupled to the
+      qubits under study (resonator_spectroscopy)
+    - Having calibrated qubits' pi and pi/2 pulses by running qubit
+      spectroscopy, rabi_chevron, power_rabi and updating the config
+    - Having calibrated the readout (readout_frequency, amplitude,
+      duration_optimization IQ_blobs) for each qubit
+      for better SNR, and having saved the derived readout threshold values in
+      the config
+    - Set the desired flux biases in the case of flux-tunable qubits
+
+This script implements the logic outlined in the notebook
+https://github.com/bornman-nick/quantum-state-and-process-tomography. See that
+notebook for a detailed explanation of standard two-qubit state tomography
+"""
+
+import numpy as np
+from scipy.linalg import sqrtm
+
+from qm.qua import *
+from qm import QuantumMachinesManager
+from qm import SimulationConfig
+from configuration import *
+from qualang_tools.results import progress_counter, fetching_tool
+
+from helper_functions import rotated_multiplexed_state_discrimination
+
+
+###################
+# The QUA program #
+###################
+
+# qubits under test, assuming there are multiple qubits
+# with XY line elements "q<qubit>_xy", flux lines
+# "q<qubit>_z", and readout resonator elements "rr<qubit>"
+
+qubit1 = 1
+qubit2 = 2
+
+if qubit1 == 1:
+    threshold1 = ge_threshold_q1
+elif qubit1 == 2:
+    threshold1 = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit1 number chosen")
+
+if qubit2 == 1:
+    threshold2 = ge_threshold_q1
+elif qubit2 == 2:
+    threshold2 = ge_threshold_q2
+else:
+    raise ValueError(f"Incorrect qubit2 number chosen")
+
+if qubit1 == qubit2:
+    raise ValueError(f"The value of qubit1 cannot equal that of qubit2")
+
+
+n_avg = 10_000
+
+
+# subroutine to prepare desired qubit state
+def prepare_state(qubit1, qubit2):
+    # write whatever QUA code you need in order to create the
+    # state to perform tomography on, from an initial
+    # ground state. For example, to create the |1> 1/sqrt(2)(|0>+i|1>)
+    # state
+    play("y180", f"q{qubit1}_xy")
+    play("-x90", f"q{qubit2}_xy")
+
+
+# matrix representation of the ideal composite qubit state from above
+# (|1> 1/sqrt(2)(|0>+i|1>)) (<1| 1/sqrt(2)(<0|-i<1|))
+ideal_state = np.array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 1 / 2, 1j / 2], [0, 0, -1j / 2, 1 / 2]])
+
+
+with program() as two_qubit_state_tomography:
+
+    n = declare(int)
+    n_st = declare_stream()
+    I = [declare(fixed) for _ in range(2)]
+    Q = [declare(fixed) for _ in range(2)]
+    # I_st = [declare_stream() for _ in range(2)]
+    # Q_st = [declare_stream() for _ in range(2)]
+
+    states = [
+        declare(bool),
+        declare(bool),
+    ]  # QUA variable for the measured qubit states
+    # states_st = [declare_stream(), declare_stream()]  # Stream for the qubits states
+
+    p00 = declare(int)  # variable to track number of |0>|0> measurements
+    p01 = declare(int)  # variable to track number of |0>|1> measurements
+    p10 = declare(int)  # variable to track number of |1>|0> measurements
+    p11 = declare(int)  # variable to track number of |1>|1> measurements
+
+    prob_vec_results_st = declare_stream()  # Stream for the probability vector - average of states vector
+
+    c = declare(int)  # QUA variable for switching between projections
+
+    with for_(n, 0, n < n_avg, n + 1):  # QUA for_ loop for averaging
+        with for_(c, 1, c <= 15, c + 1):  # QUA for_ loop for switching between projections
+
+            prepare_state(qubit1, qubit2)
+            align()
+
+            assign(p00, 0)
+            assign(p01, 0)
+            assign(p10, 0)
+            assign(p11, 0)
+
+            with switch_(c):
+                with case_(1):  # projection along Z1, X2
+
+                    play("-y90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P10 - P01 - P11
+
+                with case_(2):  # projection along Z1, Y2
+
+                    play("x90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P10 - P01 - P11
+
+                with case_(3):  # projection along Z1, Z2
+
+                    wait(pi_len // 4, f"q{qubit1}_xy", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P10 - P01 - P11
+
+                with case_(4):  # projection along X1, Z2
+
+                    play("-y90", f"q{qubit1}_xy")
+                    # Stokes parameter for this case: P00 + P01 - P10 - P11
+
+                with case_(5):  # projection along X1, X2
+
+                    play("-y90", f"q{qubit1}_xy")
+                    play("-y90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(6):  # projection along X1, Y2
+
+                    play("-y90", f"q{qubit1}_xy")
+                    play("x90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(7):  # projection along X1, Z2
+
+                    play("-y90", f"q{qubit1}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(8):  # projection along Y1, Z2
+
+                    play("x90", f"q{qubit1}_xy")
+                    # Stokes parameter for this case: P00 + P01 - P10 - P11
+
+                with case_(9):  # projection along Y1, X2
+
+                    play("x90", f"q{qubit1}_xy")
+                    play("-y90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(10):  # projection along Y1, Y2
+
+                    play("x90", f"q{qubit1}_xy")
+                    play("x90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(11):  # projection along Y1, Z2
+
+                    play("x90", f"q{qubit1}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(12):  # projection along Z1, Z2
+
+                    wait(pi_len // 4, f"q{qubit1}_xy", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P01 - P10 - P11
+
+                with case_(13):  # projection along Z1, X2
+
+                    play("-y90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(14):  # projection along Z1, Y2
+
+                    play("x90", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                with case_(15):  # projection along Z1, Z2
+
+                    wait(pi_len // 4, f"q{qubit1}_xy", f"q{qubit2}_xy")
+                    # Stokes parameter for this case: P00 + P11 - P10 - P01
+
+                align()
+
+            rotated_multiplexed_state_discrimination(
+                I,
+                None,
+                Q,
+                None,
+                states,
+                None,
+                [qubit1, qubit2],
+                [threshold1, threshold2],
+            )
+
+            align()
+
+            # populate correct variable depending on which qubit states
+            # were measured
+            with if_(states[0]):
+                with if_(states[1]):
+                    assign(p11, 1)  # |1>|1> was measured
+                with else_():
+                    assign(p10, 1)  # |1>|0> was measured
+            with else_():
+                with if_(states[1]):
+                    assign(p01, 1)  # |0>|1> was measured
+                with else_():
+                    assign(p00, 1)  # |0>|0> was measured
+
+            # stream all four possible values - only a single one of these
+            # four should be 1, the rest, 0
+            save(p00, prob_vec_results_st)
+            save(p01, prob_vec_results_st)
+            save(p10, prob_vec_results_st)
+            save(p11, prob_vec_results_st)
+
+            wait(thermalization_time * u.ns)
+
+        save(n, n_st)
+
+    with stream_processing():
+        n_st.save("iteration")
+        prob_vec_results_st.buffer(4).buffer(15).average().save("probs")
+
+
+#####################################
+#  Open Communication with the QOP  #
+#####################################
+qmm = QuantumMachinesManager(host=qop_ip, port=qop_port, cluster_name=cluster_name, octave=octave_config)
+
+###########################
+# Run or Simulate Program #
+###########################
+simulate = False
+
+if simulate:
+    # Simulates the QUA program for the specified duration
+    simulation_config = SimulationConfig(duration=10_000)  # In clock cycles = 4ns
+    job = qmm.simulate(config, two_qubit_state_tomography, simulation_config)
+    job.get_simulated_samples().con1.plot()
+
+else:
+    # Open the quantum machine
+    qm = qmm.open_qm(config)
+    # Send the QUA program to the OPX, which compiles and executes it
+    job = qm.execute(two_qubit_state_tomography)
+    # Get results from QUA program
+    results = fetching_tool(job, data_list=["probs", "iteration"], mode="live")
+
+    while results.is_processing():
+
+        # Fetch results
+        iteration, probs = results.fetch_all()
+
+        # Progress bar
+        progress_counter(iteration, n_avg, start_time=results.get_start_time())
+
+    # Close the quantum machines at the end in order to put all flux biases to
+    # 0 so that the fridge doesn't heat-up
+    qm.close()
+
+    # Use probs to reconstruct density matrix
+
+    # initialise vector to contain the 15 stokes parameters
+    stokes = np.zeros(15)
+
+    # For cases 1, 2,and 3 - Stokes parameter is P00 + P10 - P01 - P11
+    # For cases 4, 8 and 12 - Stokes parameter is P00 + P01 - P10 - P11
+    # For cases 5, 6, 7, 9, 10, 11, 13, 14, 15 - P00 + P11 - P10 - P01
+    for i in range(1, 16):
+        if i in [1, 2, 3]:
+            stokes[i - 1] = probs[i - 1][0] + probs[i - 1][2] - probs[i - 1][1] - probs[i - 1][3]
+        elif i in [4, 8, 12]:
+            stokes[i - 1] = probs[i - 1][0] + probs[i - 1][1] - probs[i - 1][2] - probs[i - 1][3]
+        elif i in [5, 6, 7, 9, 10, 11, 13, 14, 15]:
+            stokes[i - 1] = probs[i - 1][0] + probs[i - 1][3] - probs[i - 1][1] - probs[i - 1][2]
+
+    # Derive the density matrix
+    I = np.array([[1, 0], [0, 1]])
+    sigma_x = np.array([[0, 1], [1, 0]])
+    sigma_y = np.array([[0, -1j], [1j, 0]])
+    sigma_z = np.array([[1, 0], [0, -1]])
+
+    # Density matrix Pauli operator basis
+    II = np.kron(I, I)
+    IX = np.kron(I, sigma_x)
+    IY = np.kron(I, sigma_y)
+    IZ = np.kron(I, sigma_z)
+    XI = np.kron(sigma_x, I)
+    XX = np.kron(sigma_x, sigma_x)
+    XY = np.kron(sigma_x, sigma_y)
+    XZ = np.kron(sigma_x, sigma_z)
+    YI = np.kron(sigma_y, I)
+    YX = np.kron(sigma_y, sigma_x)
+    YY = np.kron(sigma_y, sigma_y)
+    YZ = np.kron(sigma_y, sigma_z)
+    ZI = np.kron(sigma_z, I)
+    ZX = np.kron(sigma_z, sigma_x)
+    ZY = np.kron(sigma_z, sigma_y)
+    ZZ = np.kron(sigma_z, sigma_z)
+
+    rho = 0.25 * (
+        II
+        + stokes[0] * IX
+        + stokes[1] * IY
+        + stokes[2] * IZ
+        + stokes[3] * XI
+        + stokes[4] * XX
+        + stokes[5] * XY
+        + stokes[6] * XZ
+        + stokes[7] * YI
+        + stokes[8] * YX
+        + stokes[9] * YY
+        + stokes[10] * YZ
+        + stokes[11] * ZI
+        + stokes[12] * ZX
+        + stokes[13] * ZY
+        + stokes[14] * ZZ
+    )
+
+    print(f"The density matrix is:\n{np.round(rho, decimals=3)}")
+
+    sqrt_ideal_state = sqrtm(ideal_state)
+
+    state_fidelity = (np.abs(sqrtm(sqrt_ideal_state @ rho @ sqrt_ideal_state).trace())) ** 2
+
+    print(f"The state fidelity is: {np.round(state_fidelity, decimals=4)}")
diff --git a/Tutorials/intro-to-simulation/cloud-simulator-example.py b/Tutorials/intro-to-simulation/cloud-simulator-example.py
new file mode 100644
index 000000000..6ad046fcf
--- /dev/null
+++ b/Tutorials/intro-to-simulation/cloud-simulator-example.py
@@ -0,0 +1,73 @@
+from qm import QuantumMachinesManager, SimulationConfig
+from qm.qua import play, program
+from qm_saas import QoPVersion, QmSaas
+
+# for qm-saas version 1.1.1
+
+# Define quantum machine configuration dictionary
+config = {
+    "version": 1,
+    "controllers": {
+        "con1": {
+            "analog_outputs": {
+                1: {"offset": +0.0},
+            },
+        }
+    },
+    "elements": {
+        "qe1": {
+            "singleInput": {"port": ("con1", 1)},
+            "intermediate_frequency": 5e6,
+            "operations": {
+                "playOp": "constPulse",
+            },
+        },
+    },
+    "pulses": {
+        "constPulse": {
+            "operation": "control",
+            "length": 1000,  # in ns
+            "waveforms": {"single": "const_wf"},
+        },
+    },
+    "waveforms": {
+        "const_wf": {"type": "constant", "sample": 0.2},
+    },
+}
+
+# ======================================================================================================================
+# Description: Use default host and port for the QOP simulator. Email and password are mandatory.
+#              Using the context manager ensures them simulator instance is closed properly.
+# ======================================================================================================================
+
+# These should be changed to your credentials.
+email = "john.doe@mail.com"
+password = "Password_given_by_QM"
+
+# Initialize QOP simulator client
+client = QmSaas(email=email, password=password)
+
+# Choose your QOP version (QOP2.x.y or QOP3.x.y)
+# if you choose QOP3.x.y, make sure you are using an adequate config.
+version = QoPVersion.v2_2_2
+
+with client.simulator(version=version) as instance:  # Specify the QOP version
+    # Initialize QuantumMachinesManager with the simulation instance details
+    qmm = QuantumMachinesManager(
+        host=instance.host, port=instance.port, connection_headers=instance.default_connection_headers
+    )
+
+    # Define a QUA program
+    with program() as prog:
+        play("playOp", "qe1")
+
+    # Open quantum machine with the provided configuration and simulate the QUA program
+    qm = qmm.open_qm(config)
+    job = qm.simulate(prog, SimulationConfig(int(1000)))
+
+    # Retrieve and handle simulated samples
+    samples = job.get_simulated_samples()
+    print("Test passed")
+
+# analysis of simulation
+# do something with samples
diff --git a/Tutorials/intro-to-simulation/readme.md b/Tutorials/intro-to-simulation/readme.md
index f8ca6d3a4..c9f000055 100644
--- a/Tutorials/intro-to-simulation/readme.md
+++ b/Tutorials/intro-to-simulation/readme.md
@@ -5,41 +5,57 @@ slug: ./
 id: index
 ---
 
-This example shows usage of the hardware simulator. The simulator mimics the output of the hardware, 
-with its exact timing and voltage level, following the compilation of QUA to the FPGA's low-level. 
-It is a useful tool for predicting and debugging the outcome of a QUA program. 
-Read more on the simulator, and it's capabilities in the [QUA docs](https://qm-docs.qualang.io/guides/simulator).
-
-The Examples
+This script intro-to-simulation.py shows usage of the hardware simulator. The simulator mimics the output of the
+hardware,
+with its exact timing and voltage level, following the compilation of QUA to the FPGA's low-level.
+It is a useful tool for predicting and debugging the outcome of a QUA program.
+Read more on the simulator, and it's capabilities in
+the [QUA docs](https://docs.quantum-machines.co/latest/docs/Guides/simulator/).
+
+In addition, the script cloud-simulator-example.py demonstrates the usage
+of [QM cloud simulator](https://docs.quantum-machines.co/latest/docs/Guides/qm_saas_guide/), enabling to run the HW
+simulator
+on virtual instances (and doesn't require access to actual HW).
+
+The intro-to-simulation Examples
 ============
 
-The following script gives three examples for the usage of the simulator in the QOP.
+The script gives three examples for the usage of the simulator in the QOP.
 
-First and second exampled
+First and second examples
 -------------
 
-The first and second demonstrates basic usage od the simulator. Note that the simulation is done via the QuantumMachinesManager (qmm) object. 
-The simulation function is given the configuration, the program and the `SimulationConfig` instance. The last sets the duration for the 
-simulation. i.e. how many clock cycles should the simulator simulate. The outcome of the simulator is a `job` object. 
+The first and second examples demonstrate basic usage of the simulator. Note that the simulation is done via the
+QuantumMachinesManager (qmm) object.
+The simulation function is given the configuration, the program and the `SimulationConfig` instance. The last sets the
+duration for the
+simulation. i.e. how many clock cycles should the simulator simulate. The outcome of the simulator is a `job` object.
 
-The examples also demonstrate how to obtain and plot the simulated output of the hardware. Finally, the first example also demonstrates
-how the simulator can simulate saving variables to a stream, as would occur in the real hardware. 
+The examples also demonstrate how to obtain and plot the simulated output of the hardware. Finally, the first example
+also demonstrates
+how the simulator can simulate saving variables to a stream, as would occur in the real hardware.
 
 Third Example
 ==============
 
 The third example demonstrates a slightly more advanced usage. It shows how a loop-back connection can be defined,
-to simulate acquisition and demodulation of ADC signals. In the example, a connection from analog output 1 of controller 1
-is connected to analog input 1 of controller 1. The example also shows that the demodulation and adc input can be simulated
-and saved to the stream processing, which later can be fetched and analyzed. It is important to note that the data will only
-be available if the simulation duration was long enough to simulate it. In the current example, to "fill" the buffer in the stream processing,
-the simulation must simulate the entire program with all the loop's iterations.  
+to simulate acquisition and demodulation of ADC signals. In the example, a connection from analog output 1 of controller
+1
+is connected to analog input 1 of controller 1. The example also shows that the demodulation and adc input can be
+simulated
+and saved to the stream processing, which later can be fetched and analyzed. It is important to note that the data will
+only
+be available if the simulation duration was long enough to simulate it. In the current example, to "fill" the buffer in
+the stream processing,
+the simulation must simulate the entire program with all the loop's iterations.
 
 Fourth Example
 =============
 
-The fourth example demonstrates a simulation of a multi-controllers system. This is done by specifying the connectivity between the different controllers.
-This is important since The exact timing of multi-controllers operations is dependent on that connectivity configuration.
+The fourth example demonstrates a simulation of a multi-controllers system. This is done by specifying the connectivity
+between the different controllers.
+This is important since The exact timing of multi-controllers operations is dependent on that connectivity
+configuration.
 In the example we use a tool to create the controllers' connections in the format that is required by the simulator.
 The tool is available in our (very useful) repo [py-qua-tools](https://github.com/qua-platform/py-qua-tools).
 
diff --git a/pyproject.toml b/pyproject.toml
index 11897f7c2..4b7a2eb01 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,6 +17,7 @@ flake8 = "^3.9.1"
 
 [tool.black]
 line-length = 120
+#force-exclude = '''Two-Flux-Tunable-Transmons/Use Case 4 - Single- and Two-Qubit State and Process Tomography/helper_functions.py'''
 #include = '(examples)/.+\.pyi?$'
 #include = '(Quantum Control Applications)/.+\.pyi?$'
 #include = '(Tutorials)/.+\.pyi?$'