BTD: remove old/legacy back-transformed diagnostics (#3485)

* Start removing old BTD

* Remove GetCellCenteredData

* Remove do_backtransform_fields and do_backtransform_particles

* Remove more functions

* Remove more variables

* Update documentation

* Fix CI test `RigidInjection_BTD`

* Remove slicing from `BTD_ReducedSliceDiag`

* Rename `BTD_ReducedSliceDiag` as `LaserAcceleration_BTD`

* Query deprecated input and abort

Co-authored-by: Edoardo Zoni <ezoni@lbl.gov>

29 files changed, 85 insertions, 2671 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 84add7daa..21f7e7db4 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -955,11 +955,6 @@ Particle initialization
         lead to memory issues if not periodically cleared. To clear the buffer
         call ``warpx_clearParticleBoundaryBuffer()``.
 
-* ``<species>.do_back_transformed_diagnostics`` (`0` or `1` optional, default `1`)
-    Only used when ``warpx.do_back_transformed_diagnostics=1``. When running in a
-    boosted frame, whether or not to plot back-transformed diagnostics for
-    this species.
-
 * ``<species>.do_field_ionization`` (`0` or `1`) optional (default `0`)
     Do field ionization for this species (using the ADK theory).
 
@@ -2182,81 +2177,6 @@ BackTransformed Diagnostics (with support for Plotfile/openPMD output)
     value for buffer size and use slices to reduce the memory footprint and maintain
     optimum I/O performance.
 
-Back-Transformed Diagnostics (legacy output)
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-
-``BackTransformedDiagnostics`` are used when running a simulation in a boosted frame, to reconstruct output data to the lab frame
-
-* ``warpx.do_back_transformed_diagnostics`` (`0` or `1`)
-    Whether to use the **back-transformed diagnostics** (i.e. diagnostics that
-    perform on-the-fly conversion to the laboratory frame, when running
-    boosted-frame simulations)
-
-* ``warpx.lab_data_directory`` (`string`)
-    The directory in which to save the lab frame data when using the
-    **back-transformed diagnostics**. If not specified, the default is
-    is `lab_frame_data`.
-
-* ``warpx.num_snapshots_lab`` (`integer`)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1``.
-    The number of lab-frame snapshots that will be written.
-
-* ``warpx.dt_snapshots_lab`` (`float`, in seconds)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1``.
-    The time interval inbetween the lab-frame snapshots (where this
-    time interval is expressed in the laboratory frame).
-
-* ``warpx.dz_snapshots_lab`` (`float`, in meters)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1``.
-    Distance between the lab-frame snapshots (expressed in the laboratory
-    frame). ``dt_snapshots_lab`` is then computed by
-    ``dt_snapshots_lab = dz_snapshots_lab/c``. Either `dt_snapshots_lab`
-    or `dz_snapshot_lab` is required.
-
-* ``warpx.do_back_transformed_fields`` (`0 or 1`)
-    Whether to use the **back-transformed diagnostics** for the fields.
-
-* ``warpx.back_transformed_diag_fields`` (space-separated list of `string`)
-    Which fields to dumped in back-transformed diagnostics. Choices are
-    'Ex', 'Ey', Ez', 'Bx', 'By', Bz', 'jx', 'jy', jz' and 'rho'. Example:
-    ``warpx.back_transformed_diag_fields = Ex Ez By``. By default, all fields
-    are dumped.
-
-* ``warpx.buffer_size`` (`integer`)
-    The default size of the back transformed diagnostic buffers used to generate lab-frame
-    data is 256. That is, when the multifab with lab-frame data has 256 z-slices,
-    the data will be flushed out. However, if many lab-frame snapshots are required for
-    diagnostics and visualization, the GPU may run out of memory with many large boxes with
-    a size of 256 in the z-direction. This input parameter can then be used to set a
-    smaller buffer-size, preferably multiples of 8, such that, a large number of
-    lab-frame snapshot data can be generated without running out of gpu memory.
-    The downside to using a small buffer size, is that the I/O time may increase due
-    to frequent flushes of the lab-frame data. The other option is to keep the default
-    value for buffer size and use slices to reduce the memory footprint and maintain
-    optimum I/O performance.
-
-* ``slice.num_slice_snapshots_lab`` (`integer`)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1``.
-    The number of back-transformed field and particle data that
-    will be written for the reduced domain defined by ``slice.dom_lo``
-    and ``slice.dom_hi``. Note that the 'slice' is a reduced
-    diagnostic which could be 1D, 2D, or 3D, aligned with the co-ordinate axes.
-    These slices can be visualized using read_raw_data.py and the HDF5 format can
-    be visualized using the h5py library. Please see the documentation on visualization
-    for further details.
-
-* ``slice.dt_slice_snapshots_lab`` (`float`, in seconds)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1``.
-    The time interval between the back-transformed reduced diagnostics (where this
-    time interval is expressed in the laboratory frame).
-
-* ``slice.particle_slice_width_lab`` (`float`, in meters)
-    Only used when ``warpx.do_back_transformed_diagnostics`` is ``1`` and
-    ``slice.num_slice_snapshots_lab`` is non-zero. Particles are
-    copied from the full back-transformed diagnostic to the reduced
-    slice diagnostic if there are within the user-defined width from
-    the slice region defined by ``slice.dom_lo`` and ``slice.dom_hi``.
-
 Boundary Scraping Diagnostics
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py
index 8881eac7b..ccb551832 100755
--- a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py
+++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py
@@ -25,7 +25,6 @@ import sys
 
 import numpy as np
 import openpmd_api as io
-import read_raw_data
 from scipy.constants import m_e
 import yt
 
@@ -36,18 +35,10 @@ import checksumAPI
 
 filename = sys.argv[1]
 
-# Tolerances to check consistency between legacy BTD and new BTD
+# Tolerances to check consistency between plotfile BTD and openPMD BTD
 rtol = 1e-16
 atol = 1e-16
 
-# Read data from legacy back-transformed diagnostics
-snapshot = './lab_frame_data/snapshots/snapshot00001'
-x_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'x')
-z_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'z')
-ux_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'ux')
-uy_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'uy')
-uz_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'uz')
-
 # Read data from new back-transformed diagnostics (plotfile)
 ds_plotfile = yt.load(filename)
 x_plotfile = ds_plotfile.all_data()['beam', 'particle_position_x'].v
@@ -66,19 +57,12 @@ uy_openpmd = ds_openpmd.particles['beam']['momentum']['y'][:]
 uz_openpmd = ds_openpmd.particles['beam']['momentum']['z'][:]
 series.flush()
 
-# Sort and compare arrays to check consistency between legacy BTD and new BTD (plotfile)
-assert(np.allclose(np.sort(x_legacy), np.sort(x_plotfile), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(z_legacy), np.sort(z_plotfile), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(ux_legacy*m_e), np.sort(ux_plotfile), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(uy_legacy*m_e), np.sort(uy_plotfile), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(uz_legacy*m_e), np.sort(uz_plotfile), rtol=rtol, atol=atol))
-
-# Sort and compare arrays to check consistency between legacy BTD and new BTD (openPMD)
-assert(np.allclose(np.sort(x_legacy), np.sort(x_openpmd), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(z_legacy), np.sort(z_openpmd), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(ux_legacy*m_e), np.sort(ux_openpmd), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(uy_legacy*m_e), np.sort(uy_openpmd), rtol=rtol, atol=atol))
-assert(np.allclose(np.sort(uz_legacy*m_e), np.sort(uz_openpmd), rtol=rtol, atol=atol))
+# Sort and compare arrays to check consistency between plotfile BTD and openPMD BTD
+assert(np.allclose(np.sort(x_plotfile), np.sort(x_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(z_plotfile), np.sort(z_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(ux_plotfile), np.sort(ux_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uy_plotfile), np.sort(uy_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uz_plotfile), np.sort(uz_openpmd), rtol=rtol, atol=atol))
 
 # Initial parameters
 z0 = 20.e-6
@@ -86,8 +70,8 @@ x0 = 1.e-6
 theta0 = np.arcsin(0.1)
 
 # Theoretical beam width after propagation with rigid injection
-z = np.mean(z_legacy)
-x = np.std(x_legacy)
+z = np.mean(z_plotfile)
+x = np.std(x_plotfile)
 print(f'Beam position = {z}')
 print(f'Beam width    = {x}')
 
diff --git a/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame b/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame
index d5c4c2b3c..9bf858b8d 100644
--- a/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame
+++ b/Examples/Modules/RigidInjection/inputs_2d_BoostedFrame
@@ -68,8 +68,3 @@ diag2.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
 diag2.format = openpmd
 diag2.openpmd_backend = h5
 diag2.buffer_size = 32
-
-# old BTD diagnostics
-warpx.do_back_transformed_diagnostics = 1
-warpx.num_snapshots_lab = 2
-warpx.dt_snapshots_lab = 1.8679589331096515e-13
diff --git a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py b/Examples/Modules/boosted_diags/analysis.py
index 6fa4e9c93..c6c089f98 100755
--- a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py
+++ b/Examples/Modules/boosted_diags/analysis.py
@@ -21,7 +21,6 @@ import sys
 
 import numpy as np
 import openpmd_api as io
-import read_raw_data
 import yt
 
 yt.funcs.mylog.setLevel(0)
@@ -35,22 +34,6 @@ filename = sys.argv[1]
 rtol = 1e-16
 atol = 1e-16
 
-# Read data from legacy back-transformed diagnostics (entire domain)
-snapshot = './lab_frame_data/snapshots/snapshot00003'
-header   = './lab_frame_data/snapshots/Header'
-allrd, info = read_raw_data.read_lab_snapshot(snapshot, header)
-Ez_legacy = allrd['Ez']
-print(f'Ez_legacy.shape = {Ez_legacy.shape}')
-Ez_legacy_1D = np.squeeze(Ez_legacy[Ez_legacy.shape[0]//2,Ez_legacy.shape[1]//2,:])
-
-# Read data from reduced back-transformed diagnostics (slice)
-snapshot_slice = './lab_frame_data/slices/slice00003'
-header_slice   = './lab_frame_data/slices/Header'
-allrd, info = read_raw_data.read_lab_snapshot(snapshot_slice, header_slice)
-Ez_legacy_slice = allrd['Ez']
-print(f'Ez_legacy_slice.shape = {Ez_legacy_slice.shape}')
-Ez_legacy_slice_1D = np.squeeze(Ez_legacy_slice[Ez_legacy_slice.shape[0]//2,1,:])
-
 # Read data from new back-transformed diagnostics (plotfile)
 ds_plotfile = yt.load(filename)
 data = ds_plotfile.covering_grid(
@@ -69,12 +52,5 @@ series.flush()
 # Compare arrays to check consistency between new BTD formats (plotfile and openPMD)
 assert(np.allclose(Ez_plotfile, Ez_openpmd, rtol=rtol, atol=atol))
 
-# Check slicing
-err = np.max(np.abs(Ez_legacy_slice_1D-Ez_legacy_1D)) / np.max(np.abs(Ez_legacy_1D))
-tol = 1e-16
-print(f'error = {err}')
-print(f'tolerance = {tol}')
-assert(err < tol)
-
 test_name = os.path.split(os.getcwd())[1]
 checksumAPI.evaluate_checksum(test_name, filename)
diff --git a/Examples/Modules/boosted_diags/inputs_3d_slice b/Examples/Modules/boosted_diags/inputs_3d
index ff68b5909..ba98558be 100644
--- a/Examples/Modules/boosted_diags/inputs_3d_slice
+++ b/Examples/Modules/boosted_diags/inputs_3d
@@ -55,7 +55,6 @@ electrons.zmax = .003
 electrons.profile = constant
 electrons.density = 3.5e24
 electrons.do_continuous_injection = 1
-electrons.do_back_transformed_diagnostics = 1
 
 ions.charge = q_e
 ions.mass = m_p
@@ -71,7 +70,6 @@ ions.zmax = .003
 ions.profile = constant
 ions.density = 3.5e24
 ions.do_continuous_injection = 1
-ions.do_back_transformed_diagnostics = 1
 
 beam.charge = -q_e
 beam.mass = m_e
@@ -104,14 +102,6 @@ laser1.profile_t_peak = 40.e-15    # The time at which the laser reaches its pea
 laser1.profile_focal_distance = 0.5e-3  # Focal distance from the antenna (in meters)
 laser1.wavelength = 0.81e-6         # The wavelength of the laser (in meters)
 
-slice.dom_lo     = xmin  0.0  zmin
-slice.dom_hi     = xmax  0.0  zmax
-slice.coarsening_ratio = 1 1 1
-slice.plot_int = -1
-slice.num_slice_snapshots_lab = 4
-slice.dt_slice_snapshots_lab = 3.3356409519815207e-12
-slice.particle_slice_width_lab = 2.e-6
-
 # Diagnostics
 diagnostics.diags_names = diag1 diag2
 
@@ -132,9 +122,3 @@ diag2.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
 diag2.format = openpmd
 diag2.buffer_size = 32
 diag2.openpmd_backend = h5
-
-# old BTD diagnostics
-warpx.do_back_transformed_diagnostics = 1
-warpx.num_snapshots_lab = 4
-warpx.dz_snapshots_lab = 0.001
-warpx.back_transformed_diag_fields= Ex Ey Ez By rho
diff --git a/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost b/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
index b2c76e685..76dcd3ee2 100644
--- a/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
+++ b/Examples/Physics_applications/plasma_acceleration/inputs_2d_boost
@@ -38,9 +38,6 @@ algo.particle_shape = 3
 #################################
 warpx.gamma_boost = 10.0
 warpx.boost_direction = z
-warpx.do_back_transformed_diagnostics = 1
-warpx.num_snapshots_lab = 22
-warpx.dt_snapshots_lab = 3.335640951981521e-11
 
 #################################
 ############ PLASMA #############
diff --git a/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost b/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
index c1c630dca..2264872ec 100644
--- a/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
+++ b/Examples/Physics_applications/plasma_acceleration/inputs_3d_boost
@@ -38,9 +38,6 @@ algo.particle_shape = 3
 #################################
 warpx.gamma_boost = 10.0
 warpx.boost_direction = z
-warpx.do_back_transformed_diagnostics = 1
-warpx.num_snapshots_lab = 22
-warpx.dt_snapshots_lab = 3.335640951981521e-11
 
 #################################
 ############ PLASMA #############
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index c112f0157..30b9f4c7c 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -1922,9 +1922,6 @@ class LabFrameFieldDiagnostic(picmistandard.PICMI_LabFrameFieldDiagnostic,
 
     Parameters
     ----------
-    warpx_new_BTD: bool, optional
-        Use the new BTD diagnostics
-
     warpx_format: string, optional
         Passed to <diagnostic name>.format
 
@@ -1947,36 +1944,16 @@ class LabFrameFieldDiagnostic(picmistandard.PICMI_LabFrameFieldDiagnostic,
         Passed to <diagnostic name>.upper_bound
     """
     def init(self, kw):
-        self.use_new_BTD = kw.pop('warpx_new_BTD', False)
-        if self.use_new_BTD:
-            # The user is using the new BTD
-            self.format = kw.pop('warpx_format', None)
-            self.openpmd_backend = kw.pop('warpx_openpmd_backend', None)
-            self.file_prefix = kw.pop('warpx_file_prefix', None)
-            self.file_min_digits = kw.pop('warpx_file_min_digits', None)
-            self.buffer_size = kw.pop('warpx_buffer_size', None)
-            self.lower_bound = kw.pop('warpx_lower_bound', None)
-            self.upper_bound = kw.pop('warpx_upper_bound', None)
+        # The user is using the new BTD
+        self.format = kw.pop('warpx_format', None)
+        self.openpmd_backend = kw.pop('warpx_openpmd_backend', None)
+        self.file_prefix = kw.pop('warpx_file_prefix', None)
+        self.file_min_digits = kw.pop('warpx_file_min_digits', None)
+        self.buffer_size = kw.pop('warpx_buffer_size', None)
+        self.lower_bound = kw.pop('warpx_lower_bound', None)
+        self.upper_bound = kw.pop('warpx_upper_bound', None)
 
     def initialize_inputs(self):
-        if self.use_new_BTD:
-            self.initialize_inputs_new()
-        else:
-            self.initialize_inputs_old()
-
-    def initialize_inputs_old(self):
-
-        pywarpx.warpx.check_consistency('num_snapshots_lab', self.num_snapshots, 'The number of snapshots must be the same in all lab frame diagnostics')
-        pywarpx.warpx.check_consistency('dt_snapshots_lab', self.dt_snapshots, 'The time between snapshots must be the same in all lab frame diagnostics')
-        pywarpx.warpx.check_consistency('lab_data_directory', self.write_dir, 'The write directory must be the same in all lab frame diagnostics')
-
-        pywarpx.warpx.do_back_transformed_diagnostics = 1
-        pywarpx.warpx.num_snapshots_lab = self.num_snapshots
-        pywarpx.warpx.dt_snapshots_lab = self.dt_snapshots
-        pywarpx.warpx.do_back_transformed_fields = 1
-        pywarpx.warpx.lab_data_directory = self.write_dir
-
-    def initialize_inputs_new(self):
 
         self.add_diagnostic()
 
@@ -2025,30 +2002,6 @@ class LabFrameFieldDiagnostic(picmistandard.PICMI_LabFrameFieldDiagnostic,
 
         self.set_write_dir()
 
-
-class LabFrameParticleDiagnostic(picmistandard.PICMI_LabFrameParticleDiagnostic):
-    def initialize_inputs(self):
-
-        pywarpx.warpx.check_consistency('num_snapshots_lab', self.num_snapshots, 'The number of snapshots must be the same in all lab frame diagnostics')
-        pywarpx.warpx.check_consistency('dt_snapshots_lab', self.dt_snapshots, 'The time between snapshots must be the same in all lab frame diagnostics')
-        pywarpx.warpx.check_consistency('lab_data_directory', self.write_dir, 'The write directory must be the same in all lab frame diagnostics')
-
-        pywarpx.warpx.do_back_transformed_diagnostics = 1
-
-        if isinstance(self.species, Species):
-            self.species.do_back_transformed_diagnostics = 1
-        else:
-            try:
-                for specie in self.species:
-                    specie.do_back_transformed_diagnostics = 1
-            except TypeError:
-                pass
-
-        pywarpx.warpx.num_snapshots_lab = self.num_snapshots
-        pywarpx.warpx.dt_snapshots_lab = self.dt_snapshots
-        pywarpx.warpx.lab_data_directory = self.write_dir
-
-
 class ReducedDiagnostic(picmistandard.base._ClassWithInit, WarpXDiagnosticBase):
     """
     Sets up a reduced diagnostic in the simulation.
diff --git a/Regression/Checksum/benchmarks_json/BTD_ReducedSliceDiag.json b/Regression/Checksum/benchmarks_json/LaserAcceleration_BTD.json
index 1fd9e23a3..1fd9e23a3 100644
--- a/Regression/Checksum/benchmarks_json/BTD_ReducedSliceDiag.json
+++ b/Regression/Checksum/benchmarks_json/LaserAcceleration_BTD.json
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 94fe785d7..6f2f038cb 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -231,12 +231,11 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 doComparison = 0
-aux1File = Tools/PostProcessing/read_raw_data.py
 analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py
 
-[BTD_ReducedSliceDiag]
+[LaserAcceleration_BTD]
 buildDir = .
-inputFile = Examples/Modules/boosted_diags/inputs_3d_slice
+inputFile = Examples/Modules/boosted_diags/inputs_3d
 runtime_params =
 dim = 3
 addToCompileString = USE_OPENPMD=TRUE
@@ -250,8 +249,7 @@ compileTest = 0
 doVis = 0
 compareParticles = 0
 doComparison = 0
-aux1File = Tools/PostProcessing/read_raw_data.py
-analysisRoutine = Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py
+analysisRoutine = Examples/Modules/boosted_diags/analysis.py
 
 [nci_corrector]
 buildDir = .
diff --git a/Source/Diagnostics/BackTransformedDiagnostic.H b/Source/Diagnostics/BackTransformedDiagnostic.H
deleted file mode 100644
index deb32e7a6..000000000
--- a/Source/Diagnostics/BackTransformedDiagnostic.H
+++ /dev/null
@@ -1,287 +0,0 @@
-/* Copyright 2019 Andrew Myers, Axel Huebl, Maxence Thevenet
- * Revathi Jambunathan, Weiqun Zhang
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#ifndef WARPX_BackTransformedDiagnostic_H_
-#define WARPX_BackTransformedDiagnostic_H_
-
-#include "BackTransformedDiagnostic_fwd.H"
-
-#include "Particles/MultiParticleContainer.H"
-#include "Particles/WarpXParticleContainer.H"
-
-#include <AMReX_Box.H>
-#include <AMReX_IntVect.H>
-#include <AMReX_MultiFab.H>
-#include <AMReX_REAL.H>
-#include <AMReX_RealBox.H>
-#include <AMReX_Vector.H>
-
-#include <AMReX_BaseFwd.H>
-
-#include <map>
-#include <memory>
-#include <string>
-#include <vector>
-
-/** \brief
- * The capability for back-transformed lab-frame data is implemented to generate
- * the full diagnostic snapshot for the entire domain and reduced diagnostic
- * (1D, 2D or 3D 'slices') for a sub-domain.
- * LabFrameDiag class defines the parameters required to backtrasform data from
- * boosted frame at (z_boost,t_boost) to lab-frame at (z_lab, t_lab) using Lorentz
- * transformation. This Lorentz transformation picks out one slice corresponding
- * to both of those times, at position current_z_boost and current_z_lab in the
- * boosted and lab frames, respectively.
- * Two derived classes, namely, LabFrameSnapShot and LabFrameSlice are defined to
- * store the full back-transformed diagnostic snapshot of the entire domain and
- * reduced back-transformed diagnostic for a sub-domain, respectively.
- * The approach here is to define an array of LabFrameDiag  which would include
- * both, full domain snapshots and reduced domain 'slices', sorted based on their
- * respective t_lab. This is done to re-use the backtransformed data stored in
- * the slice multifab at (z_lab,t_lab)
- * for the full domain snapshot and sub-domain slices that have the same t_lab,
- * instead of re-generating the backtransformed slice data at z_lab for a given
- * t_lab for each diagnostic.
- */
-class LabFrameDiag {
-    public:
-    std::string m_file_name;
-    amrex::Real m_t_lab;
-    amrex::RealBox m_prob_domain_lab_;
-    amrex::IntVect m_prob_ncells_lab_;
-    amrex::RealBox m_diag_domain_lab_;
-    amrex::Box m_buff_box_;
-
-    amrex::Real m_current_z_lab;
-    amrex::Real m_current_z_boost;
-    amrex::Real m_inv_gamma_boost_;
-    amrex::Real m_inv_beta_boost_;
-    amrex::Real m_dz_lab_;
-    amrex::Real m_particle_slice_dx_lab_;
-
-    int m_ncomp_to_dump_;
-    std::vector<std::string> m_mesh_field_names_;
-
-    int m_file_num;
-
-    // For back-transformed diagnostics of grid fields, data_buffer_
-    // stores a buffer of the fields in the lab frame (in a MultiFab, i.e.
-    // with all box data etc.). When the buffer if full, dump to file.
-    std::unique_ptr<amrex::MultiFab> m_data_buffer_;
-    // particles_buffer_ is currently blind to refinement level.
-    // particles_buffer_[j] is a WarpXParticleContainer::DiagnosticParticleData
-    //  where - j is the species number for the current diag
-    amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> m_particles_buffer_;
-    // buff_counter_ is the number of z slices in data_buffer_
-    int m_buff_counter_;
-    int m_num_buffer_ = 256;
-    int m_max_box_size = 256;
-    void updateCurrentZPositions(amrex::Real t_boost, amrex::Real inv_gamma,
-                                 amrex::Real inv_beta);
-
-    void createLabFrameDirectories();
-
-    void writeLabFrameHeader();
-
-    /// Back-transformed lab-frame field data is copied from
-    /// tmp_slice to data_buffer where it is stored.
-    /// For the full diagnostic, all the data in the MultiFab is copied.
-    /// For the reduced diagnostic, the data is selectively copied if the
-    /// extent of the z_lab multifab intersects with the user-defined sub-domain
-    /// for the reduced diagnostic (i.e., a 1D, 2D, or 3D region of the domain).
-    virtual void AddDataToBuffer(amrex::MultiFab& /*tmp_slice_ptr*/, int /*i_lab*/,
-                 amrex::Vector<int> const& /*map_actual_fields_to_dump*/){}
-
-    /// Back-transformed lab-frame particles is copied from
-    /// tmp_particle_buffer to particles_buffer.
-    /// For the full diagnostic, all the particles are copied,
-    /// while for the reduced diagnostic, particles are selectively
-    /// copied if their position in within the user-defined
-    /// sub-domain +/- 1 cell size width for the reduced slice diagnostic.
-    virtual void AddPartDataToParticleBuffer(
-                 amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> const& /*tmp_particle_buffer*/,
-                 int /*nSpeciesBoostedFrame*/) {}
-
-    // The destructor should also be a virtual function, so that
-    // a pointer to subclass of `LabFrameDiag` actually calls the subclass's destructor.
-    virtual ~LabFrameDiag() = default;
-};
-
-/** \brief
- * LabFrameSnapShot stores the back-transformed lab-frame metadata
- * corresponding to a single time snapshot of the full domain.
- * The snapshot data is written to disk in the directory lab_frame_data/snapshots/.
- * zmin_lab, zmax_lab, and t_lab are all constant for a given snapshot.
- * current_z_lab and current_z_boost for each snapshot are updated as the
- * simulation time in the boosted frame advances.
- */
-
-class LabFrameSnapShot : public LabFrameDiag {
-    public:
-    LabFrameSnapShot(amrex::Real t_lab_in, amrex::Real t_boost,
-                     amrex::Real inv_gamma_boost_in, amrex::Real inv_beta_boost_in,
-                     amrex::Real dz_lab_in, amrex::RealBox prob_domain_lab,
-                     amrex::IntVect prob_ncells_lab, int ncomp_to_dump,
-                     std::vector<std::string> mesh_field_names,
-                     amrex::RealBox diag_domain_lab,
-                     amrex::Box diag_box, int file_num_in, const int max_box_size,
-                     const int buffer_size);
-    void AddDataToBuffer( amrex::MultiFab& tmp_slice, int k_lab,
-         amrex::Vector<int> const& map_actual_fields_to_dump) override;
-    void AddPartDataToParticleBuffer(
-         amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> const& tmp_particle_buffer,
-         int nSpeciesBoostedFrame) override;
-};
-
-
-/** \brief
- * LabFrameSlice stores the back-transformed metadata corresponding to a single time at the
- * user-defined slice location. This could be a 1D line, 2D slice, or 3D box
- * (a reduced back-transformed diagnostic) within the computational
- * domain, as defined in the input file by the user. The slice is written to disk in the
- * lab_frame_data/slices. Similar to snapshots, zmin_lab, zmax_lab, and
- * t_lab are constant for a given slice. current_z_lab and current_z_boost
- * for each snapshot are updated as the sim time in boosted frame advances.
- */
-class LabFrameSlice : public LabFrameDiag {
-    public:
-    LabFrameSlice(amrex::Real t_lab_in, amrex::Real t_boost,
-                     amrex::Real inv_gamma_boost_in, amrex::Real inv_beta_boost_in,
-                     amrex::Real dz_lab_in, amrex::RealBox prob_domain_lab,
-                     amrex::IntVect prob_ncells_lab, int ncomp_to_dump,
-                     std::vector<std::string> mesh_field_names,
-                     amrex::RealBox diag_domain_lab,
-                     amrex::Box diag_box, int file_num_in,
-                     amrex::Real particle_slice_dx_lab,
-                     const int max_box_size,
-                     const int buffer_size);
-    void AddDataToBuffer( amrex::MultiFab& tmp_slice_ptr, int i_lab,
-         amrex::Vector<int> const& map_actual_fields_to_dump) override;
-    void AddPartDataToParticleBuffer(
-         amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> const& tmp_particle_buffer,
-         int nSpeciesBoostedFrame) override;
-};
-
-/** \brief
- * BackTransformedDiagnostic class handles the back-transformation of data when
- * running simulations in a boosted frame of reference to the lab-frame.
- * Because of the relativity of simultaneity, events that are synchronized
- * in the simulation boosted frame are not
- * synchronized in the lab frame. Thus, at a given t_boost, we must write
- * slices of back-transformed data to multiple output files, each one
- * corresponding to a given time in the lab frame. The member function
- * writeLabFrameData() orchestrates the operations required to
- * Lorentz-transform data from boosted-frame to lab-frame and store them
- * in the LabFrameDiag class, which writes out the field and particle data
- * to the output directory. The functions Flush() and writeLabFrameData()
- * are called at the end of the simulation and when the
- * the buffer for data storage is full, respectively. The particle data
- * is collected and written only if particle.do_back_transformed_diagnostics = 1.
- */
-class BackTransformedDiagnostic {
-
-public:
-
-    BackTransformedDiagnostic (amrex::Real zmin_lab, amrex::Real zmax_lab,
-                               amrex::Real v_window_lab, amrex::Real dt_snapshots_lab,
-                               int N_snapshots, amrex::Real dt_slice_snapshots_lab,
-                               int N_slice_snapshots, amrex::Real gamma_boost,
-                               amrex::Real t_boost, amrex::Real dt_boost,
-                               int boost_direction, const amrex::Geometry& geom,
-                               amrex::RealBox& slice_realbox,
-                               amrex::Real particle_slice_width_lab);
-
-    /// Flush() is called at the end of the simulation when the buffers that contain
-    /// back-transformed lab-frame data even if they are not full.
-    void Flush (const amrex::Geometry& geom);
-
-    /** \brief
-     * The order of operations performed in writeLabFrameData is as follows :
-     * 1. Loops over the sorted back-transformed diags and for each diag
-     *    steps 2-7 are performed
-     * 2. Based on t_lab and t_boost, obtain z_lab and z_boost.
-     * 3. Define data_buffer multifab that will store the data in the BT diag.
-     * 4. Define slice multifab at z_index that corresponds to z_boost and
-     *    getslicedata using cell-centered data at z_index and its distribution map.
-     * 5. Lorentz transform data stored in slice from z_boost,t_Boost to z_lab,t_lab
-     *    and store in slice multifab.
-     * 6. Generate a temporary slice multifab with distribution map of lab-frame
-     *    data but at z_boost and
-     *    ParallelCopy data from the slice multifab to the temporary slice.
-     * 7. Finally, AddDataToBuffer is called where the data from temporary slice
-     *    is simply copied from tmp_slice(i,j,k_boost) to
-     *    LabFrameDiagSnapshot(i,j,k_lab) for full BT lab-frame diagnostic
-     *    OR from tmp_slice(i,j,k_boost) to
-     *    LabFrameDiagSlice(i,j,k_lab) for the reduced slice diagnostic
-     * 8. Similarly, particles that crossed the z_boost plane are selected
-     *    and lorentz-transformed to the lab-frame and copied to the full
-     *    and reduce diagnostic and stored in particle_buffer.
-     */
-    void writeLabFrameData (const amrex::MultiFab* cell_centered_data,
-                            const MultiParticleContainer& mypc,
-                            const amrex::Geometry& geom,
-                            const amrex::Real t_boost, const amrex::Real dt);
-    /// The metadata containg information on t_boost, num_snapshots, and Lorentz parameters.
-    void writeMetaData ();
-
-private:
-    amrex::Real m_gamma_boost_;
-    amrex::Real m_inv_gamma_boost_;
-    amrex::Real m_beta_boost_;
-    amrex::Real m_inv_beta_boost_;
-    amrex::Real m_dz_lab_;
-    amrex::Real m_inv_dz_lab_;
-    amrex::Real m_dt_snapshots_lab_;
-    amrex::Real m_dt_boost_;
-    int m_N_snapshots_;
-    int m_boost_direction_;
-    int m_N_slice_snapshots_;
-    amrex::Real m_dt_slice_snapshots_lab_;
-    amrex::Real m_particle_slice_width_lab_;
-
-    int m_num_buffer_ = 256;
-    int m_max_box_size_ = 256;
-
-    std::vector<std::unique_ptr<LabFrameDiag> > m_LabFrameDiags_;
-
-    void writeParticleData (
-         const WarpXParticleContainer::DiagnosticParticleData& pdata,
-         const std::string& name, const int i_lab);
-
-#ifdef WARPX_USE_HDF5
-    void writeParticleDataHDF5(
-         const WarpXParticleContainer::DiagnosticParticleData& pdata,
-         const std::string& name, const std::string& species_name);
-#endif
-    // Map field names and component number in cell_centered_data
-    std::map<std::string, int> m_possible_fields_to_dump = {
-        {"Ex" , 0},
-        {"Ey" , 1},
-        {"Ez" , 2},
-        {"Bx" , 3},
-        {"By" , 4},
-        {"Bz" , 5},
-        {"jx" , 6},
-        {"jy" , 7},
-        {"jz" , 8},
-        {"rho", 9} };
-
-    // maps field index in data_buffer_[i] -> cell_centered_data for
-    // snapshots i. By default, all fields in cell_centered_data are dumped.
-    // Needs to be amrex::Vector because used in a ParallelFor kernel.
-    amrex::Vector<int> map_actual_fields_to_dump;
-    // Name of fields to dump. By default, all fields in cell_centered_data.
-    // Needed for file headers only.
-    std::vector<std::string> m_mesh_field_names = {"Ex", "Ey", "Ez",
-                                                   "Bx", "By", "Bz",
-                                                   "jx", "jy", "jz", "rho"};
-    int m_ncomp_to_dump = 10;
-
-
-};
-
-#endif
diff --git a/Source/Diagnostics/BackTransformedDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp
deleted file mode 100644
index 7d148abf1..000000000
--- a/Source/Diagnostics/BackTransformedDiagnostic.cpp
+++ /dev/null
@@ -1,1662 +0,0 @@
-/* Copyright 2019 Andrew Myers, Axel Huebl, Maxence Thevenet
- * Revathi Jambunathan, Weiqun Zhang
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-#include "BackTransformedDiagnostic.H"
-
-#include "Utils/Parser/ParserUtils.H"
-#include "Utils/TextMsg.H"
-#include "Utils/WarpXConst.H"
-#include "Utils/WarpXProfilerWrapper.H"
-#include "Utils/TextMsg.H"
-#include "WarpX.H"
-
-#include <ablastr/utils/Communication.H>
-
-#include <AMReX_Array4.H>
-#include <AMReX_BLassert.H>
-#include <AMReX_BoxArray.H>
-#include <AMReX_Config.H>
-#include <AMReX_DistributionMapping.H>
-#include <AMReX_Extension.H>
-#include <AMReX_FArrayBox.H>
-#include <AMReX_FabArray.H>
-#include <AMReX_Geometry.H>
-#include <AMReX_GpuContainers.H>
-#include <AMReX_GpuControl.H>
-#include <AMReX_GpuLaunch.H>
-#include <AMReX_GpuQualifiers.H>
-#include <AMReX_MFIter.H>
-#include <AMReX_MultiFabUtil.H>
-#include <AMReX_PODVector.H>
-#include <AMReX_ParallelDescriptor.H>
-#include <AMReX_ParmParse.H>
-#include <AMReX_PlotFileUtil.H>
-#include <AMReX_SPACE.H>
-#include <AMReX_Scan.H>
-#include <AMReX_StructOfArrays.H>
-#include <AMReX_Utility.H>
-#include <AMReX_VectorIO.H>
-#include <AMReX_VisMF.H>
-
-#ifdef WARPX_USE_HDF5
-    #include <hdf5.h>
-#endif
-
-#include <algorithm>
-#include <cmath>
-#include <fstream>
-#include <memory>
-
-using namespace amrex;
-
-namespace
-{
-    constexpr int permission_flag_rwxrxrx = 0755;
-}
-
-#ifdef WARPX_USE_HDF5
-
-/*
-  Helper functions for doing the HDF5 IO.
-
- */
-namespace
-{
-
-    const std::vector<std::string> particle_field_names = {"w", "x", "y",
-                                                           "z", "ux", "uy", "uz"};
-
-    /*
-      Creates the HDF5 file in truncate mode and closes it.
-      Should be run only by the root process.
-    */
-    void output_create (const std::string& file_path) {
-        WARPX_PROFILE("output_create");
-        hid_t file = H5Fcreate(file_path.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            file >=0,
-            "Error: could not create file at " + file_path
-        )
-        H5Fclose(file);
-    }
-
-    /*
-      Writes a single string attribute to the given group.
-      Should only be called by the root process.
-    */
-    void write_string_attribute (hid_t& group, const std::string& key, const std::string& val)
-    {
-        hid_t str_type     = H5Tcopy(H5T_C_S1);
-        hid_t scalar_space = H5Screate(H5S_SCALAR);
-
-        // Fix the str_type length for the format string.
-        H5Tset_size(str_type, strlen(val.c_str()));
-
-        hid_t attr = H5Acreate(group, key.c_str(), str_type, scalar_space, H5P_DEFAULT, H5P_DEFAULT);
-        H5Awrite(attr, str_type, val.c_str());
-
-        H5Aclose(attr);
-        H5Sclose(scalar_space);
-        H5Tclose(str_type);
-    }
-
-    /*
-      Writes a single double attribute to the given group.
-      Should only be called by the root process.
-    */
-    void write_double_attribute (hid_t& group, const std::string& key, const double val)
-    {
-        hid_t scalar_space = H5Screate(H5S_SCALAR);
-
-        hid_t attr = H5Acreate(group, key.c_str(), H5T_IEEE_F32LE, scalar_space,
-                               H5P_DEFAULT, H5P_DEFAULT);
-        H5Awrite(attr, H5T_NATIVE_DOUBLE, &val);
-
-        H5Aclose(attr);
-        H5Sclose(scalar_space);
-    }
-
-    /*
-      Opens the output file and writes all of metadata attributes.
-      Should be run only by the root process.
-    */
-    void output_write_metadata (const std::string& file_path,
-                                const int istep, const Real time, const Real dt)
-    {
-        WARPX_PROFILE("output_write_metadata");
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-
-        write_string_attribute(file, "software", "warpx");
-        write_string_attribute(file, "softwareVersion", "0.0.0");
-        write_string_attribute(file, "meshesPath", "fields/");
-        write_string_attribute(file, "iterationEncoding", "fileBased");
-        write_string_attribute(file, "iterationFormat", "data%T.h5");
-        write_string_attribute(file, "openPMD", "1.1.0");
-        write_string_attribute(file, "basePath", "/data/%T/");
-
-        hid_t group = H5Gcreate(file, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        group = H5Gcreate(group, std::to_string(istep).c_str(),
-                          H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-        write_double_attribute(group, "time", time);
-        write_double_attribute(group, "timeUnitSI", 1.0);
-        write_double_attribute(group, "dt", dt);
-
-        // Field groups
-        group = H5Gcreate(group, "fields", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-        // Close all resources.
-        H5Gclose(group);
-        H5Fclose(file);
-        H5close();
-    }
-
-    /*
-      Creates a dataset with the given cell dimensions, at the path
-      "/native_fields/(field_name)".
-      Should be run only by the root rank.
-    */
-    void output_create_field (const std::string& file_path, const std::string& field_path,
-                              const unsigned nx, const unsigned ny, const unsigned nz)
-    {
-        WARPX_PROFILE("output_create_field");
-
-        // Open the output.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-        // Create a 3D, nx x ny x nz dataspace.
-#if defined(WARPX_DIM_3D)
-        hsize_t dims[3] = {nx, ny, nz};
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-        hsize_t dims[3] = {nx, nz};
-#else
-        hsize_t dims[3] = {nz};
-#endif
-        hid_t grid_space = H5Screate_simple(AMREX_SPACEDIM, dims, NULL);
-
-        // Create the dataset.
-        hid_t dataset = H5Dcreate(file, field_path.c_str(), H5T_IEEE_F64LE,
-                                  grid_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            dataset >=0,
-            "Error: could not create dataset. H5 returned "
-            + std::to_string(dataset))
-        );
-
-        // Close resources.
-        H5Dclose(dataset);
-        H5Sclose(grid_space);
-        H5Fclose(file);
-    }
-
-    /*
-      Creates a group associated with a single particle species.
-      Should be run by all processes collectively.
-    */
-    void output_create_species_group (const std::string& file_path, const std::string& species_name)
-    {
-        MPI_Comm comm = MPI_COMM_WORLD;
-        MPI_Info info = MPI_INFO_NULL;
-        int mpi_rank;
-        MPI_Comm_rank(comm, &mpi_rank);
-
-        // Create the file access prop list.
-        hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS);
-        H5Pset_fapl_mpio(pa_plist, comm, info);
-
-        // Open the output.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist);
-
-        hid_t group = H5Gcreate(file, species_name.c_str(),
-                                H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-        H5Gclose(group);
-        H5Fclose(file);
-
-    }
-
-    /*
-      Resize an extendible dataset, suitable for storing particle data.
-      Should be run only by the root rank.
-    */
-    long output_resize_particle_field (const std::string& file_path, const std::string& field_path,
-                                       const long num_to_add)
-    {
-        WARPX_PROFILE("output_resize_particle_field");
-
-        // Open the output.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-
-        int rank;
-        hsize_t dims[1];
-
-        hid_t dataset = H5Dopen2 (file, field_path.c_str(), H5P_DEFAULT);
-        hid_t filespace = H5Dget_space (dataset);
-        rank = H5Sget_simple_extent_ndims (filespace);
-        herr_t status = H5Sget_simple_extent_dims (filespace, dims, NULL);
-
-        // set new size
-        hsize_t new_size[1];
-        new_size[0] = dims[0] + num_to_add;
-        status = H5Dset_extent (dataset, new_size);
-
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            status >= 0,
-            "Error: set extent filed on dataset "
-            + std::to_string(dataset)
-        );
-
-        // Close resources.
-        H5Sclose(filespace);
-        H5Dclose(dataset);
-        H5Fclose(file);
-
-        return dims[0];
-    }
-
-    /*
-      Writes to a dataset that has been extended to the proper size. Suitable for writing particle data.
-      Should be run on all ranks collectively.
-    */
-    void output_write_particle_field (const std::string& file_path, const std::string& field_path,
-                                      const Real* data_ptr, const long count, const long index)
-    {
-        WARPX_PROFILE("output_write_particle_field");
-
-        MPI_Comm comm = MPI_COMM_WORLD;
-        MPI_Info info = MPI_INFO_NULL;
-        int mpi_rank;
-        MPI_Comm_rank(comm, &mpi_rank);
-
-        // Create the file access prop list.
-        hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS);
-        H5Pset_fapl_mpio(pa_plist, comm, info);
-
-        // Open the output.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist);
-
-        int RANK = 1;
-        hsize_t offset[1];
-        hsize_t dims[1];
-        herr_t status;
-
-        hid_t dataset = H5Dopen (file, field_path.c_str(), H5P_DEFAULT);
-
-        // Make sure the dataset is there.
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            dataset >= 0,
-            "Error on rank " + std::to_string(mpi_rank)
-            +". Count not find dataset " + field_path
-        );
-
-        hid_t filespace = H5Dget_space (dataset);
-
-        offset[0] = index;
-        dims[0] = count;
-
-        // Create collective io prop list.
-        hid_t collective_plist = H5Pcreate(H5P_DATASET_XFER);
-        H5Pset_dxpl_mpio(collective_plist, H5FD_MPIO_INDEPENDENT);
-
-        if (count > 0) {
-
-            /* Define memory space */
-            hid_t memspace = H5Screate_simple (RANK, dims, NULL);
-
-            status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL,
-                                          dims, NULL);
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                status >= 0,
-                "Error on rank " + std::to_string(ParallelDescriptor::MyProc())
-                +" could not select hyperslab."
-            );
-
-            /* Write the data to the extended portion of dataset  */
-            status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace,
-                              filespace, collective_plist, data_ptr);
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                status >= 0,
-                "Error on rank " + std::to_string(ParallelDescriptor::MyProc())
-                +" could not write hyperslab."
-            );
-
-            status = H5Sclose (memspace);
-        }
-
-        ParallelDescriptor::Barrier();
-
-        // Close resources.
-        H5Pclose(collective_plist);
-        H5Sclose(filespace);
-        H5Dclose(dataset);
-        H5Fclose(file);
-        H5Pclose(pa_plist);
-    }
-
-    /*
-      Creates an extendible dataset, suitable for storing particle data.
-      Should be run on all ranks collectively.
-    */
-    void output_create_particle_field (const std::string& file_path, const std::string& field_path)
-    {
-        WARPX_PROFILE("output_create_particle_field");
-
-        MPI_Comm comm = MPI_COMM_WORLD;
-        MPI_Info info = MPI_INFO_NULL;
-        int mpi_rank;
-        MPI_Comm_rank(comm, &mpi_rank);
-
-        // Create the file access prop list.
-        hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS);
-        H5Pset_fapl_mpio(pa_plist, comm, info);
-
-        // Open the output.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist);
-
-        constexpr int RANK = 1;
-        hsize_t dims[1] = {0};
-        hsize_t maxdims[1] = {H5S_UNLIMITED};
-        hsize_t      chunk_dims[2] = {4};
-
-        hid_t dataspace = H5Screate_simple (RANK, dims, maxdims);
-
-        // Enable chunking
-        hid_t prop = H5Pcreate (H5P_DATASET_CREATE);
-        herr_t status = H5Pset_chunk (prop, RANK, chunk_dims);
-
-        hid_t dataset = H5Dcreate2 (file, field_path.c_str(), H5T_NATIVE_DOUBLE, dataspace,
-                                    H5P_DEFAULT, prop, H5P_DEFAULT);
-
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            dataset >= 0,
-            "Error: could not create dataset. H5 returned "
-            + std::to_string(dataset)
-        );
-
-        // Close resources.
-        H5Dclose(dataset);
-        H5Pclose(prop);
-        H5Sclose(dataspace);
-        H5Fclose(file);
-    }
-
-    /*
-      Write the only component in the multifab to the dataset given by field_name.
-      Uses hdf5-parallel.
-    */
-    void output_write_field (const std::string& file_path,
-                             const std::string& field_path,
-                             const MultiFab& mf, const int comp,
-                             const int lo_x, const int lo_y, const int lo_z)
-    {
-
-        WARPX_PROFILE("output_write_field");
-
-        MPI_Comm comm = MPI_COMM_WORLD;
-        MPI_Info info = MPI_INFO_NULL;
-        int mpi_rank;
-        MPI_Comm_rank(comm, &mpi_rank);
-
-        // Create the file access prop list.
-        hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS);
-        H5Pset_fapl_mpio(pa_plist, comm, info);
-
-        // Open the file, and the group.
-        hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist);
-        // Open the field dataset.
-        hid_t dataset = H5Dopen(file, field_path.c_str(), H5P_DEFAULT);
-
-        // Make sure the dataset is there.
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            dataset >= 0,
-            "Error on rank " + std::to_string(mpi_rank)
-            +". Count not find dataset " + field_path
-        );
-
-        // Grab the dataspace of the field dataset from file.
-        hid_t file_dataspace = H5Dget_space(dataset);
-
-        // Create collective io prop list.
-        hid_t collective_plist = H5Pcreate(H5P_DATASET_XFER);
-        H5Pset_dxpl_mpio(collective_plist, H5FD_MPIO_INDEPENDENT);
-
-        // Iterate over Fabs, select matching hyperslab and write.
-        hid_t status;
-        // slab lo index and shape.
-#if defined(WARPX_DIM_3D)
-        hsize_t slab_offsets[3], slab_dims[3];
-        int shift[3];
-        shift[0] = lo_x;
-        shift[1] = lo_y;
-        shift[2] = lo_z;
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-        hsize_t slab_offsets[2], slab_dims[2];
-        int shift[2];
-        shift[0] = lo_x;
-        shift[1] = lo_z;
-#else
-        hsize_t slab_offsets[1], slab_dims[1];
-        int shift[1];
-        shift[0] = lo_z;
-#endif
-        hid_t slab_dataspace;
-
-        int write_count = 0;
-
-        std::vector<Real> transposed_data;
-
-        for (MFIter mfi(mf); mfi.isValid(); ++mfi)
-        {
-            const Box& box = mfi.validbox();
-            const int *lo_vec = box.loVect();
-            const int *hi_vec = box.hiVect();
-
-            transposed_data.resize(box.numPts(), 0.0);
-
-            // Set slab offset and shape.
-            for (int idim = 0; idim < AMREX_SPACEDIM; ++idim)
-            {
-                AMREX_ASSERT(lo_vec[idim] >= 0);
-                AMREX_ASSERT(hi_vec[idim] > lo_vec[idim]);
-                slab_offsets[idim] = lo_vec[idim] - shift[idim];
-                slab_dims[idim] = hi_vec[idim] - lo_vec[idim] + 1;
-            }
-
-            int cnt = 0;
-            AMREX_D_TERM(
-                         for (int i = lo_vec[0]; i <= hi_vec[0]; ++i),
-                         for (int j = lo_vec[1]; j <= hi_vec[1]; ++j),
-                         for (int k = lo_vec[2]; k <= hi_vec[2]; ++k))
-                transposed_data[cnt++] = mf[mfi](IntVect(AMREX_D_DECL(i, j, k)), comp);
-
-            // Create the slab space.
-            slab_dataspace = H5Screate_simple(AMREX_SPACEDIM, slab_dims, NULL);
-
-            // Select the hyperslab matching this fab.
-            status = H5Sselect_hyperslab(file_dataspace, H5S_SELECT_SET,
-                                         slab_offsets, NULL, slab_dims, NULL);
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                status >= 0,
-                "Error on rank " + std::to_string(mpi_rank)
-                +" could not select hyperslab.\n"
-            );
-
-            // Write this pencil.
-            status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, slab_dataspace,
-                              file_dataspace, collective_plist, transposed_data.data());
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                status >= 0,
-                "Error on rank " + std::to_string(mpi_rank)
-                +" could not write hyperslab."
-            );
-
-            H5Sclose(slab_dataspace);
-            write_count++;
-        }
-
-        ParallelDescriptor::Barrier();
-
-        // Close HDF5 resources.
-        H5Pclose(collective_plist);
-        H5Sclose(file_dataspace);
-        H5Dclose(dataset);
-        H5Fclose(file);
-        H5Pclose(pa_plist);
-    }
-}
-#endif
-
-bool compare_tlab_uptr (const std::unique_ptr<LabFrameDiag>&a,
-                        const std::unique_ptr<LabFrameDiag>&b)
-{
-    return a->m_t_lab < b->m_t_lab;
-}
-
-namespace
-{
-void
-LorentzTransformZ (MultiFab& data, Real gamma_boost, Real beta_boost)
-{
-    // Loop over tiles/boxes and in-place convert each slice from boosted
-    // frame to back-transformed lab frame.
-#ifdef AMREX_USE_OMP
-#pragma omp parallel if (Gpu::notInLaunchRegion())
-#endif
-    for (MFIter mfi(data, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-        const Box& tile_box = mfi.tilebox();
-        Array4< Real > arr = data[mfi].array();
-        // arr(x,y,z,comp) where 0->9 comps are
-        // Ex Ey Ez Bx By Bz jx jy jz rho
-        Real clight = PhysConst::c;
-        ParallelFor(tile_box,
-            [=] AMREX_GPU_DEVICE (int i, int j, int k)
-            {
-                // Transform the transverse E and B fields. Note that ez and bz are not
-                // changed by the tranform.
-                Real e_lab = 0.0_rt, b_lab = 0.0_rt, j_lab = 0.0_rt, r_lab = 0.0_rt;
-                e_lab = gamma_boost * (arr(i, j, k, 0) +
-                                       beta_boost*clight*arr(i, j, k, 4));
-                b_lab = gamma_boost * (arr(i, j, k, 4) +
-                                       beta_boost*arr(i, j, k, 0)/clight);
-
-                arr(i, j, k, 0) = e_lab;
-                arr(i, j, k, 4) = b_lab;
-
-                e_lab = gamma_boost * (arr(i, j, k, 1) -
-                                       beta_boost*clight*arr(i, j, k, 3));
-                b_lab = gamma_boost * (arr(i, j, k, 3) -
-                                       beta_boost*arr(i, j, k, 1)/clight);
-
-                arr(i, j, k, 1) = e_lab;
-                arr(i, j, k, 3) = b_lab;
-
-                // Transform the charge and current density. Only the z component of j is affected.
-                const int j_comp_index = 8;
-                const int r_comp_index = 9;
-
-                j_lab = gamma_boost*(arr(i, j, k, j_comp_index) +
-                                     beta_boost*clight*arr(i, j, k, j_comp_index));
-                r_lab = gamma_boost*(arr(i, j, k, r_comp_index) +
-                                     beta_boost*arr(i, j, k, r_comp_index)/clight);
-
-                arr(i, j, k, j_comp_index) = j_lab;
-                arr(i, j, k, r_comp_index) = r_lab;
-            }
-        );
-    }
-}
-}
-
-BackTransformedDiagnostic::
-BackTransformedDiagnostic (Real zmin_lab, Real zmax_lab, Real v_window_lab,
-                           Real dt_snapshots_lab, int N_snapshots,
-                           Real dt_slice_snapshots_lab, int N_slice_snapshots,
-                           Real gamma_boost, Real t_boost, Real dt_boost,
-                           int boost_direction, const Geometry& geom,
-                           amrex::RealBox& slice_realbox,
-                           amrex::Real particle_slice_width_lab)
-    : m_gamma_boost_(gamma_boost),
-      m_dt_snapshots_lab_(dt_snapshots_lab),
-      m_dt_boost_(dt_boost),
-      m_N_snapshots_(N_snapshots),
-      m_boost_direction_(boost_direction),
-      m_N_slice_snapshots_(N_slice_snapshots),
-      m_dt_slice_snapshots_lab_(dt_slice_snapshots_lab),
-      m_particle_slice_width_lab_(particle_slice_width_lab)
-{
-
-#ifdef WARPX_DIM_RZ
-    amrex::Abort(Utils::TextMsg::Err("BackTransformed diagnostics is currently not supported with RZ"));
-#endif
-    WARPX_PROFILE("BackTransformedDiagnostic::BackTransformedDiagnostic");
-
-    AMREX_ALWAYS_ASSERT(WarpX::do_back_transformed_fields or
-                        WarpX::do_back_transformed_particles);
-
-    m_inv_gamma_boost_ = 1.0_rt / m_gamma_boost_;
-    m_beta_boost_ = std::sqrt(1.0_rt - m_inv_gamma_boost_*m_inv_gamma_boost_);
-    m_inv_beta_boost_ = 1.0_rt / m_beta_boost_;
-
-    m_dz_lab_ = PhysConst::c * m_dt_boost_ * m_inv_beta_boost_ * m_inv_gamma_boost_;
-    m_inv_dz_lab_ = 1.0_rt / m_dz_lab_;
-    int Nz_lab = static_cast<int>((zmax_lab - zmin_lab) * m_inv_dz_lab_);
-#if (AMREX_SPACEDIM >= 2)
-    int Nx_lab = geom.Domain().length(0);
-#endif
-#if defined(WARPX_DIM_3D)
-    int Ny_lab = geom.Domain().length(1);
-    IntVect prob_ncells_lab = {Nx_lab, Ny_lab, Nz_lab};
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-    // Ny_lab = 1;
-    IntVect prob_ncells_lab = {Nx_lab, Nz_lab};
-#else
-    // Nx_lab = 1;
-    // Ny_lab = 1;
-    IntVect prob_ncells_lab(Nz_lab);
-#endif
-    writeMetaData();
-
-    // Query fields to dump
-    std::vector<std::string> user_fields_to_dump;
-    ParmParse pp_warpx("warpx");
-    bool do_user_fields = false;
-    do_user_fields = pp_warpx.queryarr("back_transformed_diag_fields", user_fields_to_dump);
-    if (utils::parser::queryWithParser(pp_warpx, "buffer_size", m_num_buffer_)) {
-        if (m_max_box_size_ < m_num_buffer_) m_max_box_size_ = m_num_buffer_;
-    }
-    // If user specifies fields to dump, overwrite ncomp_to_dump,
-    // map_actual_fields_to_dump and mesh_field_names.
-    for (int i = 0; i < 10; ++i)
-        map_actual_fields_to_dump.push_back(i);
-
-    if (do_user_fields){
-        m_ncomp_to_dump = static_cast<int>(user_fields_to_dump.size());
-        map_actual_fields_to_dump.resize(m_ncomp_to_dump);
-        m_mesh_field_names.resize(m_ncomp_to_dump);
-        for (int i=0; i<m_ncomp_to_dump; i++){
-            std::string fieldstr = user_fields_to_dump[i];
-            m_mesh_field_names[i] = fieldstr;
-            map_actual_fields_to_dump[i] = m_possible_fields_to_dump[fieldstr];
-        }
-    }
-
-    // allocating array with total number of lab frame diags (snapshots+slices)
-    m_LabFrameDiags_.resize(N_snapshots+N_slice_snapshots);
-
-    for (int i = 0; i < N_snapshots; ++i) {
-        // steps + initial box shift
-        Real const zmax_boost = geom.ProbHi(AMREX_SPACEDIM-1);
-        Real const t_lab =
-            i * m_dt_snapshots_lab_ +
-            m_gamma_boost_ * m_beta_boost_ * zmax_boost/PhysConst::c;
-
-        // Get simulation domain physical coordinates (in boosted frame).
-        RealBox prob_domain_lab = geom.ProbDomain();
-        // Replace z bounds by lab-frame coordinates
-        // x and y bounds are the same for back-transformed lab frame and boosted frame
-        prob_domain_lab.setLo(WARPX_ZINDEX, zmin_lab + v_window_lab * t_lab);
-        prob_domain_lab.setHi(WARPX_ZINDEX, zmax_lab + v_window_lab * t_lab);
-        Box diag_box = geom.Domain();
-        m_LabFrameDiags_[i] = std::make_unique<LabFrameSnapShot>(t_lab, t_boost,
-                                m_inv_gamma_boost_, m_inv_beta_boost_, m_dz_lab_,
-                                prob_domain_lab, prob_ncells_lab,
-                                m_ncomp_to_dump, m_mesh_field_names, prob_domain_lab,
-                                diag_box, i, m_max_box_size_, m_num_buffer_);
-    }
-
-
-    for (int i = 0; i < N_slice_snapshots; ++i) {
-
-
-        // To construct LabFrameSlice(), the location of lo() and hi() of the
-        // reduced diag is computed using the user-defined values of the
-        // reduced diag (1D, 2D, or 3D). For visualization of the diagnostics,
-        // the number of cells in each dimension is required and
-        // is computed below for the reduced back-transformed lab-frame diag,
-        // similar to the full-diag.
-        const amrex::Real* current_slice_lo = slice_realbox.lo();
-        const amrex::Real* current_slice_hi = slice_realbox.hi();
-
-        const amrex::Real zmin_slice_lab = current_slice_lo[WARPX_ZINDEX] /
-                                          ( (1._rt+m_beta_boost_)*m_gamma_boost_);
-        const amrex::Real zmax_slice_lab = current_slice_hi[WARPX_ZINDEX] /
-                                          ( (1._rt+m_beta_boost_)*m_gamma_boost_);
-        auto Nz_slice_lab = static_cast<int>(
-            (zmax_slice_lab - zmin_slice_lab) * m_inv_dz_lab_);
-#if (AMREX_SPACEDIM >= 2)
-        auto Nx_slice_lab = static_cast<int>(
-            (current_slice_hi[0] - current_slice_lo[0] ) /
-            geom.CellSize(0));
-        if (Nx_slice_lab == 0 ) Nx_slice_lab = 1;
-        // if the x-dimension is reduced, increase total_cells by 1
-        // to be consistent with the number of cells created for the output.
-        if (Nx_lab != Nx_slice_lab) Nx_slice_lab++;
-#endif
-#if defined(WARPX_DIM_3D)
-        auto Ny_slice_lab = static_cast<int>(
-            (current_slice_hi[1] - current_slice_lo[1]) /
-            geom.CellSize(1));
-        if (Ny_slice_lab == 0 ) Ny_slice_lab = 1;
-        // if the y-dimension is reduced, increase total_cells by 1
-        // to be consistent with the number of cells created for the output.
-        if (Ny_lab != Ny_slice_lab) Ny_slice_lab++;
-        amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab};
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-        amrex::IntVect slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab};
-#else
-        amrex::IntVect slice_ncells_lab(Nz_slice_lab);
-#endif
-
-        IntVect slice_lo(AMREX_D_DECL(0,0,0));
-        IntVect slice_hi(AMREX_D_DECL(1,1,1));
-
-        for ( int i_dim=0; i_dim<AMREX_SPACEDIM; ++i_dim)
-        {
-           slice_lo[i_dim] = static_cast<int>(
-               (slice_realbox.lo(i_dim) - geom.ProbLo(i_dim) -
-                0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim));
-           slice_hi[i_dim] = static_cast<int>(
-               (slice_realbox.hi(i_dim) - geom.ProbLo(i_dim) -
-                0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim));
-           if (slice_lo[i_dim] == slice_hi[i_dim])
-           {
-              slice_hi[i_dim] = slice_lo[i_dim] + 1;
-           }
-        }
-        Box stmp(slice_lo,slice_hi);
-        Box slicediag_box = stmp;
-
-        // steps + initial box shift
-        Real const zmax_boost = geom.ProbHi(AMREX_SPACEDIM-1);
-        Real const t_slice_lab =
-            i * m_dt_slice_snapshots_lab_ +
-            m_gamma_boost_ * m_beta_boost_ * zmax_boost/PhysConst::c;
-
-        RealBox prob_domain_lab = geom.ProbDomain();
-        // replace z bounds by lab-frame coordinates
-        prob_domain_lab.setLo(WARPX_ZINDEX, zmin_lab + v_window_lab * t_slice_lab);
-        prob_domain_lab.setHi(WARPX_ZINDEX, zmax_lab + v_window_lab * t_slice_lab);
-        RealBox slice_dom_lab = slice_realbox;
-        // replace z bounds of slice in lab-frame coordinates
-        // note : x and y bounds are the same for lab and boosted frames
-        // initial lab slice extent //
-        slice_dom_lab.setLo(WARPX_ZINDEX, zmin_slice_lab + v_window_lab * t_slice_lab );
-        slice_dom_lab.setHi(WARPX_ZINDEX, zmax_slice_lab +
-                                         v_window_lab * t_slice_lab );
-
-        // construct labframeslice
-        m_LabFrameDiags_[i+N_snapshots] = std::make_unique<LabFrameSlice>(t_slice_lab, t_boost,
-                                m_inv_gamma_boost_, m_inv_beta_boost_, m_dz_lab_,
-                                prob_domain_lab, slice_ncells_lab,
-                                m_ncomp_to_dump, m_mesh_field_names, slice_dom_lab,
-                                slicediag_box, i, m_particle_slice_width_lab_,
-                                m_max_box_size_, m_num_buffer_);
-    }
-    // sort diags based on their respective t_lab
-    std::stable_sort(m_LabFrameDiags_.begin(), m_LabFrameDiags_.end(), compare_tlab_uptr);
-
-    AMREX_ALWAYS_ASSERT(m_max_box_size_ >= m_num_buffer_);
-}
-
-void BackTransformedDiagnostic::Flush (const Geometry& /*geom*/)
-{
-    WARPX_PROFILE("BackTransformedDiagnostic::Flush");
-
-    VisMF::Header::Version current_version = VisMF::GetHeaderVersion();
-    VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1);
-
-    const auto & mypc = WarpX::GetInstance().GetPartContainer();
-    const std::vector<std::string> species_names = mypc.GetSpeciesNames();
-
-    // Loop over BFD snapshots
-    for (auto& lf_diags : m_LabFrameDiags_) {
-
-        Real zmin_lab = lf_diags->m_prob_domain_lab_.lo(WARPX_ZINDEX);
-        auto i_lab = static_cast<int>(
-            (lf_diags->m_current_z_lab - zmin_lab) / m_dz_lab_);
-
-        if (lf_diags->m_buff_counter_ != 0) {
-            if (WarpX::do_back_transformed_fields) {
-                const BoxArray& ba = lf_diags->m_data_buffer_->boxArray();
-                const int hi = ba[0].bigEnd(m_boost_direction_);
-                const int lo = hi - lf_diags->m_buff_counter_ + 1;
-
-                //Box buff_box = geom.Domain();
-                Box buff_box = lf_diags->m_buff_box_;
-                buff_box.setSmall(m_boost_direction_, lo);
-                buff_box.setBig(m_boost_direction_, hi);
-
-                BoxArray buff_ba(buff_box);
-                buff_ba.maxSize(m_max_box_size_);
-                DistributionMapping buff_dm(buff_ba);
-
-                const int ncomp = lf_diags->m_data_buffer_->nComp();
-
-                MultiFab tmp(buff_ba, buff_dm, ncomp, 0);
-                tmp.setVal(0.0);
-
-                ablastr::utils::communication::ParallelCopy(tmp, *lf_diags->m_data_buffer_, 0, 0, ncomp,
-                                                            IntVect(AMREX_D_DECL(0, 0, 0)),
-                                                            IntVect(AMREX_D_DECL(0, 0, 0)),
-                                                            WarpX::do_single_precision_comms);
-
-#ifdef WARPX_USE_HDF5
-                for (int comp = 0; comp < ncomp; ++comp) {
-                    output_write_field(lf_diags->m_file_name,
-                                       m_mesh_field_names[comp], tmp, comp,
-                                       lbound(buff_box).x, lbound(buff_box).y,
-                                       lbound(buff_box).z);
-                }
-#else
-                std::stringstream ss;
-                ss << lf_diags->m_file_name << "/Level_0/"
-                   << Concatenate("buffer", i_lab, 5);
-                VisMF::Write(tmp, ss.str());
-#endif
-            }
-
-            if (WarpX::do_back_transformed_particles) {
-                // Loop over species to be dumped to BFD
-                for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) {
-                    // Get species name
-                    const std::string& species_name =
-                        species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)];
-#ifdef WARPX_USE_HDF5
-                    // Dump species data
-                    writeParticleDataHDF5(lf_diags->m_particles_buffer_[j],
-                                          lf_diags->m_file_name,
-                                          species_name);
-#else
-                    std::stringstream part_ss;
-                    part_ss << lf_diags->m_file_name + "/" +
-                               species_name + "/";
-                    // Dump species data
-                    writeParticleData(lf_diags->m_particles_buffer_[j],
-                                      part_ss.str(), i_lab);
-#endif
-                }
-                lf_diags->m_particles_buffer_.clear();
-            }
-            lf_diags->m_buff_counter_ = 0;
-        }
-    }
-
-    VisMF::SetHeaderVersion(current_version);
-}
-
-
-
-
-
-void
-BackTransformedDiagnostic::
-writeLabFrameData (const MultiFab* cell_centered_data,
-                   const MultiParticleContainer& mypc,
-                   const Geometry& geom, const Real t_boost, const Real dt) {
-
-    WARPX_PROFILE("BackTransformedDiagnostic::writeLabFrameData");
-    VisMF::Header::Version current_version = VisMF::GetHeaderVersion();
-    VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1);
-
-    const RealBox& domain_z_boost = geom.ProbDomain();
-    const Real zlo_boost = domain_z_boost.lo(m_boost_direction_);
-    const Real zhi_boost = domain_z_boost.hi(m_boost_direction_);
-
-    const std::vector<std::string> species_names = mypc.GetSpeciesNames();
-    Real prev_t_lab = -dt;
-    std::unique_ptr<amrex::MultiFab> tmp_slice_ptr;
-    std::unique_ptr<amrex::MultiFab> slice;
-    amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer;
-
-    // Loop over snapshots
-    for (auto& lf_diags : m_LabFrameDiags_) {
-        // Get updated z position of snapshot
-        const Real old_z_boost = lf_diags->m_current_z_boost;
-        lf_diags->updateCurrentZPositions(t_boost,
-                                              m_inv_gamma_boost_,
-                                              m_inv_beta_boost_);
-
-        Real diag_zmin_lab = lf_diags->m_diag_domain_lab_.lo(WARPX_ZINDEX);
-        Real diag_zmax_lab = lf_diags->m_diag_domain_lab_.hi(WARPX_ZINDEX);
-
-        if ( ( lf_diags->m_current_z_boost < zlo_boost) or
-             ( lf_diags->m_current_z_boost > zhi_boost) or
-             ( lf_diags->m_current_z_lab < diag_zmin_lab) or
-             ( lf_diags->m_current_z_lab > diag_zmax_lab) ) continue;
-
-        // Get z index of data_buffer_ (i.e. in the lab frame) where
-        // simulation domain (t', [zmin',zmax']), back-transformed to lab
-        // frame, intersects with snapshot.
-        Real dom_zmin_lab = lf_diags->m_prob_domain_lab_.lo(WARPX_ZINDEX);
-        auto i_lab = static_cast<unsigned>(
-            ( lf_diags->m_current_z_lab - dom_zmin_lab) / m_dz_lab_);
-        // If buffer of snapshot i is empty...
-        if ( lf_diags->m_buff_counter_ == 0) {
-            // ... reset fields buffer data_buffer_
-            if (WarpX::do_back_transformed_fields) {
-                lf_diags->m_buff_box_.setSmall(m_boost_direction_,
-                                             i_lab - m_num_buffer_ + 1);
-                lf_diags->m_buff_box_.setBig(m_boost_direction_, i_lab);
-
-                BoxArray buff_ba(lf_diags->m_buff_box_);
-                buff_ba.maxSize(m_max_box_size_);
-                DistributionMapping buff_dm(buff_ba);
-                lf_diags->m_data_buffer_ = std::make_unique<MultiFab>(buff_ba,
-                                                buff_dm, m_ncomp_to_dump, 0);
-                lf_diags->m_data_buffer_->setVal(0.0);
-            }
-            // ... reset particle buffer particles_buffer_[i]
-            if (WarpX::do_back_transformed_particles)
-                lf_diags->m_particles_buffer_.resize(
-                                   mypc.nSpeciesBackTransformedDiagnostics());
-        }
-
-        if (WarpX::do_back_transformed_fields) {
-            const int ncomp = cell_centered_data->nComp();
-            const int start_comp = 0;
-            const bool interpolate = true;
-            // slice containing back-transformed data is generated only if t_lab != prev_t_lab and is re-used if multiple diags have the same z_lab,t_lab.
-            if (lf_diags->m_t_lab != prev_t_lab ) {
-               if (slice)
-               {
-                 slice = nullptr;
-               }
-               slice = amrex::get_slice_data(m_boost_direction_,
-                                             lf_diags->m_current_z_boost,
-                                             *cell_centered_data, geom,
-                                             start_comp, ncomp,
-                                             interpolate);
-               // Back-transform data to the lab-frame
-               LorentzTransformZ(*slice, m_gamma_boost_, m_beta_boost_);
-             }
-             // Create a 2D box for the slice in the boosted frame
-             Real dx = geom.CellSize(m_boost_direction_);
-             auto i_boost = static_cast<int>(
-                 ( lf_diags->m_current_z_boost -
-                   geom.ProbLo(m_boost_direction_))/dx);
-             //Box slice_box = geom.Domain();
-             Box slice_box = lf_diags->m_buff_box_;
-             slice_box.setSmall(m_boost_direction_, i_boost);
-             slice_box.setBig(m_boost_direction_, i_boost);
-
-             // Make it a BoxArray slice_ba
-             BoxArray slice_ba(slice_box);
-             slice_ba.maxSize(m_max_box_size_);
-             tmp_slice_ptr = std::make_unique<MultiFab>(slice_ba,
-                             lf_diags->m_data_buffer_->DistributionMap(),
-                             ncomp, 0);
-             tmp_slice_ptr->setVal(0.0);
-
-             // slice is re-used if the t_lab of a diag is equal to
-             // that of the previous diag.
-             // Back-transformed data is copied from slice
-             // which has the dmap of the domain to
-             // tmp_slice_ptr which has the dmap of the
-             // data_buffer that stores the back-transformed data.
-            ablastr::utils::communication::ParallelCopy(*tmp_slice_ptr, *slice, 0, 0, ncomp,
-                                                        IntVect(AMREX_D_DECL(0, 0, 0)),
-                                                        IntVect(AMREX_D_DECL(0, 0, 0)),
-                                                        WarpX::do_single_precision_comms);
-             lf_diags->AddDataToBuffer(*tmp_slice_ptr, i_lab,
-                                               map_actual_fields_to_dump);
-             tmp_slice_ptr = nullptr;
-        }
-
-        if (WarpX::do_back_transformed_particles) {
-
-            if (lf_diags->m_t_lab != prev_t_lab ) {
-               if (!tmp_particle_buffer.empty())
-               {
-                  tmp_particle_buffer.clear();
-                  tmp_particle_buffer.shrink_to_fit();
-               }
-               tmp_particle_buffer.resize(mypc.nSpeciesBackTransformedDiagnostics());
-               mypc.GetLabFrameData(lf_diags->m_file_name, i_lab,
-                                    m_boost_direction_, old_z_boost,
-                                    lf_diags->m_current_z_boost,
-                                    t_boost, lf_diags->m_t_lab, dt,
-                                    tmp_particle_buffer);
-            }
-            lf_diags->AddPartDataToParticleBuffer(tmp_particle_buffer,
-                               mypc.nSpeciesBackTransformedDiagnostics());
-        }
-
-        ++lf_diags->m_buff_counter_;
-        prev_t_lab = lf_diags->m_t_lab;
-        // If buffer full, write to disk.
-        if (lf_diags->m_buff_counter_ == m_num_buffer_) {
-
-            if (WarpX::do_back_transformed_fields) {
-#ifdef WARPX_USE_HDF5
-
-                Box buff_box = lf_diags->m_buff_box_;
-                for (int comp = 0; comp < lf_diags->m_data_buffer_->nComp(); ++comp)
-                    output_write_field(lf_diags->m_file_name,
-                                       m_mesh_field_names[comp],
-                                       *lf_diags->m_data_buffer_, comp,
-                                       lbound(buff_box).x, lbound(buff_box).y,
-                                       lbound(buff_box).z);
-#else
-                std::stringstream mesh_ss;
-                mesh_ss << lf_diags->m_file_name << "/Level_0/" <<
-                           Concatenate("buffer", i_lab, 5);
-                VisMF::Write( (*lf_diags->m_data_buffer_), mesh_ss.str());
-#endif
-            }
-
-            if (WarpX::do_back_transformed_particles) {
-                // Loop over species to be dumped to BFD
-                for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) {
-                    // Get species name
-                    const std::string& species_name = species_names[
-                                      mypc.mapSpeciesBackTransformedDiagnostics(j)];
-#ifdef WARPX_USE_HDF5
-                    // Write data to disk (HDF5)
-                    writeParticleDataHDF5(lf_diags->m_particles_buffer_[j],
-                                          lf_diags->m_file_name,
-                                          species_name);
-#else
-                    std::stringstream part_ss;
-
-                    part_ss << lf_diags->m_file_name + "/" +
-                               species_name + "/";
-
-                    // Write data to disk (custom)
-                    writeParticleData(lf_diags->m_particles_buffer_[j],
-                                      part_ss.str(), i_lab);
-#endif
-                }
-                lf_diags->m_particles_buffer_.clear();
-            }
-            lf_diags->m_buff_counter_ = 0;
-        }
-    }
-
-    VisMF::SetHeaderVersion(current_version);
-}
-
-#ifdef WARPX_USE_HDF5
-void
-BackTransformedDiagnostic::
-writeParticleDataHDF5 (const WarpXParticleContainer::DiagnosticParticleData& pdata,
-                       const std::string& name, const std::string& species_name)
-{
-    auto np = pdata.GetRealData(DiagIdx::w).size();
-
-    Vector<long> particle_counts(ParallelDescriptor::NProcs(), 0);
-    Vector<long> particle_offsets(ParallelDescriptor::NProcs(), 0);
-
-    ParallelAllGather::AllGather(np, particle_counts.data(),
-                                 ParallelContext::CommunicatorAll());
-
-    long total_np = 0;
-    for (int i = 0; i < ParallelDescriptor::NProcs(); ++i) {
-        particle_offsets[i] = total_np;
-        total_np += particle_counts[i];
-    }
-
-    if (total_np == 0) return;
-
-    long old_np = 0;
-    if (ParallelDescriptor::IOProcessor())
-    {
-        for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k)
-        {
-            std::string field_path = species_name + "/" + particle_field_names[k];
-            old_np = output_resize_particle_field(name, field_path, total_np);
-        }
-    }
-
-    // Note, this has the effect of an MPI Barrier between the above resize operation
-    // and the below write.
-    ParallelDescriptor::ReduceLongMax(old_np);
-
-    // Write data here
-    for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k)
-    {
-        std::string field_path = species_name + "/" + particle_field_names[k];
-        output_write_particle_field(name, field_path,
-                                    pdata.GetRealData(k).data(),
-                                    particle_counts[ParallelDescriptor::MyProc()],
-                                    particle_offsets[ParallelDescriptor::MyProc()]
-                                    + old_np);
-    }
-}
-#endif
-
-void
-BackTransformedDiagnostic::
-writeParticleData (const WarpXParticleContainer::DiagnosticParticleData& pdata,
-                   const std::string& name, const int i_lab)
-{
-    WARPX_PROFILE("BackTransformedDiagnostic::writeParticleData");
-
-    std::string field_name;
-    std::ofstream ofs;
-
-    const int MyProc = ParallelDescriptor::MyProc();
-    auto np = pdata.GetRealData(DiagIdx::w).size();
-    if (np == 0) return;
-
-    field_name = name + Concatenate("w_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::w).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("x_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::x).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("y_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::y).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("z_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::z).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("ux_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::ux).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("uy_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::uy).data(), np, ofs);
-    ofs.close();
-
-    field_name = name + Concatenate("uz_", i_lab, 5) + "_" + std::to_string(MyProc);
-    ofs.open(field_name.c_str(), std::ios::out|std::ios::binary);
-    writeData(pdata.GetRealData(DiagIdx::uz).data(), np, ofs);
-    ofs.close();
-}
-
-void
-BackTransformedDiagnostic::
-writeMetaData ()
-{
-    WARPX_PROFILE("BackTransformedDiagnostic::writeMetaData");
-
-    if (ParallelDescriptor::IOProcessor()) {
-        const std::string fullpath = WarpX::lab_data_directory + "/snapshots";
-        if (!UtilCreateDirectory(fullpath, permission_flag_rwxrxrx))
-            CreateDirectoryFailed(fullpath);
-
-        VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
-        std::ofstream HeaderFile;
-        HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
-        std::string HeaderFileName(WarpX::lab_data_directory + "/snapshots/Header");
-        HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out   |
-                                                std::ofstream::trunc |
-                                                std::ofstream::binary);
-        if(!HeaderFile.good())
-            FileOpenFailed(HeaderFileName);
-
-        HeaderFile.precision(17);
-
-        HeaderFile << m_N_snapshots_ << "\n";
-        HeaderFile << m_dt_snapshots_lab_ << "\n";
-        HeaderFile << m_gamma_boost_ << "\n";
-        HeaderFile << m_beta_boost_ << "\n";
-
-        if (m_N_slice_snapshots_ > 0) {
-           const std::string fullpath_slice = WarpX::lab_data_directory + "/slices";
-           if (!UtilCreateDirectory(fullpath_slice, permission_flag_rwxrxrx))
-               CreateDirectoryFailed(fullpath_slice);
-
-           VisMF::IO_Buffer io_buffer_slice(VisMF::IO_Buffer_Size);
-           std::ofstream HeaderFile_slice;
-           HeaderFile_slice.rdbuf()->pubsetbuf(io_buffer_slice.dataPtr(),
-                                               io_buffer_slice.size());
-           std::string HeaderFileName_slice(WarpX::lab_data_directory+
-                                            "/slices/Header");
-           HeaderFile_slice.open(HeaderFileName_slice.c_str(),
-                                                std::ofstream::out   |
-                                                std::ofstream::trunc |
-                                                std::ofstream::binary);
-
-           if (!HeaderFile_slice.good())
-               FileOpenFailed(HeaderFileName_slice);
-
-           HeaderFile_slice.precision(17);
-
-           HeaderFile_slice << m_N_slice_snapshots_ << "\n";
-           HeaderFile_slice << m_dt_slice_snapshots_lab_ << "\n";
-           HeaderFile_slice << m_gamma_boost_ << "\n";
-           HeaderFile_slice << m_beta_boost_ << "\n";
-
-        }
-
-    }
-
-
-}
-
-LabFrameSnapShot::
-LabFrameSnapShot (Real t_lab_in, Real t_boost, Real inv_gamma_boost_in,
-                  Real inv_beta_boost_in, Real dz_lab_in, RealBox prob_domain_lab,
-                  IntVect prob_ncells_lab, int ncomp_to_dump,
-                  std::vector<std::string> mesh_field_names,
-                  amrex::RealBox diag_domain_lab, Box diag_box, int file_num_in,
-                  const int max_box_size, const int num_buffer)
-{
-   m_t_lab = t_lab_in;
-   m_dz_lab_ = dz_lab_in;
-   m_inv_gamma_boost_ = inv_gamma_boost_in;
-   m_inv_beta_boost_ = inv_beta_boost_in;
-   m_prob_domain_lab_ = prob_domain_lab;
-   m_prob_ncells_lab_ = prob_ncells_lab;
-   m_diag_domain_lab_ = diag_domain_lab;
-   m_buff_box_ = diag_box;
-   m_ncomp_to_dump_ = ncomp_to_dump;
-   m_mesh_field_names_ = std::move(mesh_field_names);
-   m_file_num = file_num_in;
-   m_current_z_lab = 0.0;
-   m_current_z_boost = 0.0;
-   updateCurrentZPositions(t_boost, m_inv_gamma_boost_, m_inv_beta_boost_);
-   m_file_name = Concatenate(WarpX::lab_data_directory + "/snapshots/snapshot",
-                           m_file_num, 5);
-   createLabFrameDirectories();
-   m_buff_counter_ = 0;
-   m_max_box_size = max_box_size;
-   m_num_buffer_ = num_buffer;
-   if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr);
-}
-
-void
-LabFrameDiag::
-updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta)
-{
-    m_current_z_boost = (m_t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta;
-    m_current_z_lab   = (m_t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta;
-}
-
-void
-LabFrameDiag::
-createLabFrameDirectories() {
-#ifdef WARPX_USE_HDF5
-    if (ParallelDescriptor::IOProcessor())
-    {
-        output_create(m_file_name);
-    }
-
-    ParallelDescriptor::Barrier();
-
-    if (ParallelDescriptor::IOProcessor())
-    {
-        if (WarpX::do_back_transformed_fields)
-        {
-            const auto lo = lbound(m_buff_box_);
-            for (int comp = 0; comp < m_ncomp_to_dump_; ++comp) {
-                output_create_field(m_file_name, m_mesh_field_names_[comp],
-#if ( AMREX_SPACEDIM >= 2 )
-                                    m_prob_ncells_lab_[0],
-#else
-                                    1,
-#endif
-#if defined(WARPX_DIM_3D)
-                                    m_prob_ncells_lab_[1],
-#else
-                                    1,
-#endif
-                                    m_prob_ncells_lab_[WARPX_ZINDEX]+1);
-            }
-        }
-    }
-
-    ParallelDescriptor::Barrier();
-
-    if (WarpX::do_back_transformed_particles){
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::vector<std::string> species_names = mypc.GetSpeciesNames();
-        // Loop over species to be dumped to BFD
-        for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j)
-        {
-            // Loop over species to be dumped to BFD
-            std::string species_name =
-                species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)];
-            output_create_species_group(m_file_name, species_name);
-            for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k)
-            {
-                std::string field_path = species_name + "/" + particle_field_names[k];
-                output_create_particle_field(m_file_name, field_path);
-            }
-        }
-    }
-#else
-    if (ParallelDescriptor::IOProcessor()) {
-
-        if (!UtilCreateDirectory(m_file_name, permission_flag_rwxrxrx))
-            CreateDirectoryFailed(m_file_name);
-
-        const int nlevels = 1;
-        for(int i = 0; i < nlevels; ++i) {
-            const std::string &fullpath = LevelFullPath(i, m_file_name);
-            if (!UtilCreateDirectory(fullpath, permission_flag_rwxrxrx))
-                CreateDirectoryFailed(fullpath);
-        }
-
-        const auto & mypc = WarpX::GetInstance().GetPartContainer();
-        const std::vector<std::string> species_names = mypc.GetSpeciesNames();
-
-        const std::string particles_prefix = "particle";
-        // Loop over species to be dumped to BFD
-        for(int i = 0; i < mypc.nSpeciesBackTransformedDiagnostics(); ++i) {
-            // Get species name
-            const std::string& species_name =
-                species_names[mypc.mapSpeciesBackTransformedDiagnostics(i)];
-            const std::string fullpath = m_file_name + "/" + species_name;
-            if (!UtilCreateDirectory(fullpath, permission_flag_rwxrxrx))
-                CreateDirectoryFailed(fullpath);
-        }
-    }
-#endif
-    ParallelDescriptor::Barrier();
-
-    writeLabFrameHeader();
-}
-
-void
-LabFrameDiag::
-writeLabFrameHeader() {
-#ifndef WARPX_USE_HDF5
-    if (ParallelDescriptor::IOProcessor()) {
-        VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
-        std::ofstream HeaderFile;
-        HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
-        std::string HeaderFileName(m_file_name + "/Header");
-        HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out   |
-                                                std::ofstream::trunc |
-                                                std::ofstream::binary);
-        if(!HeaderFile.good())
-            FileOpenFailed(HeaderFileName);
-
-        HeaderFile.precision(17);
-
-        HeaderFile << m_t_lab << "\n";
-        // Write domain number of cells
-        HeaderFile << m_prob_ncells_lab_[0] << ' '
-#if defined(WARPX_DIM_3D)
-                   << m_prob_ncells_lab_[1] << ' '
-#endif
-                   << m_prob_ncells_lab_[WARPX_ZINDEX] <<'\n';
-        // Write domain physical boundaries
-        // domain lower bound
-        HeaderFile << m_diag_domain_lab_.lo(0) << ' '
-#if defined(WARPX_DIM_3D)
-                   << m_diag_domain_lab_.lo(1) << ' '
-#endif
-                   << m_diag_domain_lab_.lo(WARPX_ZINDEX) <<'\n';
-        // domain higher bound
-        HeaderFile << m_diag_domain_lab_.hi(0) << ' '
-#if defined(WARPX_DIM_3D)
-                   << m_diag_domain_lab_.hi(1) << ' '
-#endif
-                   << m_diag_domain_lab_.hi(WARPX_ZINDEX) <<'\n';
-        // List of fields dumped to file
-        for (int i=0; i<m_ncomp_to_dump_; i++)
-        {
-            HeaderFile << m_mesh_field_names_[i] << ' ';
-        }
-        HeaderFile << "\n";
-    }
-#endif
-
-}
-
-
-LabFrameSlice::
-LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in,
-                 Real inv_beta_boost_in, Real dz_lab_in, RealBox prob_domain_lab,
-                 IntVect prob_ncells_lab, int ncomp_to_dump,
-                 std::vector<std::string> mesh_field_names,
-                 RealBox diag_domain_lab, Box diag_box, int file_num_in,
-                 amrex::Real particle_slice_dx_lab, const int max_box_size,
-                 const int num_buffer)
-{
-    m_t_lab = t_lab_in;
-    m_dz_lab_ = dz_lab_in;
-    m_inv_gamma_boost_ = inv_gamma_boost_in;
-    m_inv_beta_boost_ = inv_beta_boost_in;
-    m_prob_domain_lab_ = prob_domain_lab;
-    m_prob_ncells_lab_ = prob_ncells_lab;
-    m_diag_domain_lab_ = diag_domain_lab;
-    m_buff_box_ = diag_box;
-    m_ncomp_to_dump_ = ncomp_to_dump;
-    m_mesh_field_names_ = std::move(mesh_field_names);
-    m_file_num = file_num_in;
-    m_current_z_lab = 0.0;
-    m_current_z_boost = 0.0;
-    updateCurrentZPositions(t_boost, m_inv_gamma_boost_, m_inv_beta_boost_);
-    m_file_name = Concatenate(WarpX::lab_data_directory+"/slices/slice",m_file_num,5);
-    createLabFrameDirectories();
-    m_buff_counter_ = 0;
-    m_particle_slice_dx_lab_ = particle_slice_dx_lab;
-    m_max_box_size = max_box_size;
-    m_num_buffer_ = num_buffer;
-
-    if (WarpX::do_back_transformed_fields) m_data_buffer_.reset(nullptr);
-}
-
-void
-LabFrameSnapShot::
-AddDataToBuffer( MultiFab& tmp, int k_lab,
-                 amrex::Vector<int> const& map_actual_fields_to_dump)
-{
-    const int ncomp_to_dump = static_cast<int>(map_actual_fields_to_dump.size());
-    MultiFab& buf = *m_data_buffer_;
-#ifdef AMREX_USE_GPU
-    Gpu::DeviceVector<int> d_map_actual_fields_to_dump(ncomp_to_dump);
-    Gpu::copyAsync(Gpu::hostToDevice,
-                   map_actual_fields_to_dump.begin(), map_actual_fields_to_dump.end(),
-                   d_map_actual_fields_to_dump.begin());
-    Gpu::synchronize();
-    int const* field_map_ptr = d_map_actual_fields_to_dump.dataPtr();
-#else
-    int const* field_map_ptr = map_actual_fields_to_dump.dataPtr();
-#endif
-    for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
-         Array4<Real> tmp_arr = tmp[mfi].array();
-         Array4<Real> buf_arr = buf[mfi].array();
-         // For 3D runs, tmp is a 2D (x,y) multifab that contains only
-         // slice to write to file
-         const Box& bx = mfi.tilebox();
-         ParallelFor(bx, ncomp_to_dump,
-             [=] AMREX_GPU_DEVICE(int i, int j, int k, int n)
-             {
-                 const int icomp = field_map_ptr[n];
-#if defined(WARPX_DIM_3D)
-                 buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp);
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-                 buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp);
-#else
-                 buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp);
-#endif
-             }
-         );
-    }
-}
-
-
-void
-LabFrameSlice::
-AddDataToBuffer( MultiFab& tmp, int k_lab,
-                 amrex::Vector<int> const& map_actual_fields_to_dump)
-{
-    const int ncomp_to_dump = static_cast<int>(map_actual_fields_to_dump.size());
-    MultiFab& buf = *m_data_buffer_;
-#ifdef AMREX_USE_GPU
-    Gpu::DeviceVector<int> d_map_actual_fields_to_dump(ncomp_to_dump);
-    Gpu::copyAsync(Gpu::hostToDevice,
-                   map_actual_fields_to_dump.begin(), map_actual_fields_to_dump.end(),
-                   d_map_actual_fields_to_dump.begin());
-    Gpu::synchronize();
-    int const* field_map_ptr = d_map_actual_fields_to_dump.dataPtr();
-#else
-    int const* field_map_ptr = map_actual_fields_to_dump.dataPtr();
-#endif
-    for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi)
-    {
-       const Box& bx_bf = mfi.tilebox();
-       Array4<Real> tmp_arr = tmp[mfi].array();
-       Array4<Real> buf_arr = buf[mfi].array();
-       ParallelFor(bx_bf, ncomp_to_dump,
-           [=] AMREX_GPU_DEVICE(int i, int j, int k, int n)
-           {
-              const int icomp = field_map_ptr[n];
-#if defined(WARPX_DIM_3D)
-              buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp);
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-              buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp);
-#else
-              buf_arr(k_lab,j,k,n) = tmp_arr(i,j,k,icomp);
-#endif
-           });
-    }
-
-}
-
-
-void
-LabFrameSnapShot::
-AddPartDataToParticleBuffer(
-    Vector<WarpXParticleContainer::DiagnosticParticleData> const& tmp_particle_buffer,
-    int nspeciesBoostedFrame) {
-    for (int isp = 0; isp < nspeciesBoostedFrame; ++isp) {
-        auto np = static_cast<int>(tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size());
-        if (np == 0) continue;
-
-        // allocate size of particle buffer array to np
-        // This is a growing array. Each time we add np elements
-        // to the existing array which has size = init_size
-        const int init_size = static_cast<int>(m_particles_buffer_[isp].GetRealData(DiagIdx::w).size());
-        const int total_size = init_size + np;
-        m_particles_buffer_[isp].resize(total_size);
-
-        // Data pointers to particle attributes //
-        ParticleReal* const AMREX_RESTRICT wp_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::w).data();
-        ParticleReal* const AMREX_RESTRICT x_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::x).data();
-        ParticleReal* const AMREX_RESTRICT y_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::y).data();
-        ParticleReal* const AMREX_RESTRICT z_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::z).data();
-        ParticleReal* const AMREX_RESTRICT ux_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::ux).data();
-        ParticleReal* const AMREX_RESTRICT uy_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::uy).data();
-        ParticleReal* const AMREX_RESTRICT uz_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::uz).data();
-
-        ParticleReal const* const AMREX_RESTRICT wp_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::w).data();
-        ParticleReal const* const AMREX_RESTRICT x_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::x).data();
-        ParticleReal const* const AMREX_RESTRICT y_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::y).data();
-        ParticleReal const* const AMREX_RESTRICT z_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::z).data();
-        ParticleReal const* const AMREX_RESTRICT ux_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).data();
-        ParticleReal const* const AMREX_RESTRICT uy_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).data();
-        ParticleReal const* const AMREX_RESTRICT uz_temp =
-                    tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).data();
-
-        // copy all the particles from tmp to buffer
-        amrex::ParallelFor(np,
-        [=] AMREX_GPU_DEVICE(int i)
-        {
-            wp_buff[init_size + i] = wp_temp[i];
-            x_buff[init_size + i]  = x_temp[i];
-            y_buff[init_size + i]  = y_temp[i];
-            z_buff[init_size + i]  = z_temp[i];
-            ux_buff[init_size + i] = ux_temp[i];
-            uy_buff[init_size + i] = uy_temp[i];
-            uz_buff[init_size + i] = uz_temp[i];
-        });
-    }
-}
-
-void
-LabFrameSlice::
-AddPartDataToParticleBuffer(
-    Vector<WarpXParticleContainer::DiagnosticParticleData> const& tmp_particle_buffer,
-    int nSpeciesBackTransformedDiagnostics) {
-
-
-    for (int isp = 0; isp < nSpeciesBackTransformedDiagnostics; ++isp) {
-        auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size();
-
-        if (np == 0) continue;
-
-        ParticleReal const* const AMREX_RESTRICT wp_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::w).data();
-        ParticleReal const* const AMREX_RESTRICT x_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::x).data();
-        ParticleReal const* const AMREX_RESTRICT y_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::y).data();
-        ParticleReal const* const AMREX_RESTRICT z_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::z).data();
-        ParticleReal const* const AMREX_RESTRICT ux_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).data();
-        ParticleReal const* const AMREX_RESTRICT uy_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).data();
-        ParticleReal const* const AMREX_RESTRICT uz_temp =
-             tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).data();
-
-        // temporary arrays to store copy_flag and copy_index
-        // for particles that cross the reduced domain for diagnostics.
-        amrex::Gpu::DeviceVector<int> FlagForPartCopy(np);
-        amrex::Gpu::DeviceVector<int> IndexForPartCopy(np);
-
-        int* const AMREX_RESTRICT Flag = FlagForPartCopy.dataPtr();
-        int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr();
-
-        // Compute extent of the reduced domain +/- user-defined physical width
-#if (AMREX_SPACEDIM >= 2)
-        Real const xmin = m_diag_domain_lab_.lo(0)-m_particle_slice_dx_lab_;
-        Real const xmax = m_diag_domain_lab_.hi(0)+m_particle_slice_dx_lab_;
-#endif
-#if defined(WARPX_DIM_3D)
-        Real const ymin = m_diag_domain_lab_.lo(1)-m_particle_slice_dx_lab_;
-        Real const ymax = m_diag_domain_lab_.hi(1)+m_particle_slice_dx_lab_;
-#endif
-
-        //Flag particles that need to be copied if they are
-        // within the reduced slice +/- user-defined physical width
-        amrex::ParallelFor(np,
-        [=] AMREX_GPU_DEVICE(int i)
-        {
-            Flag[i] = 0;
-#if (AMREX_SPACEDIM >= 2)
-            if ( x_temp[i] >= (xmin) &&
-                 x_temp[i] <= (xmax) )
-#endif
-            {
-#if defined(WARPX_DIM_3D)
-               if (y_temp[i] >= (ymin) &&
-                   y_temp[i] <= (ymax) )
-#endif
-               {
-                   Flag[i] = 1;
-               }
-            }
-        });
-
-        // Call exclusive scan to obtain location indices using
-        // flag values. These location indices are used to copy data
-        // from src to dst when the copy-flag is set to 1.
-        const int copy_size = amrex::Scan::ExclusiveSum(np, Flag, IndexLocation);
-        const auto init_size = static_cast<int>(
-            m_particles_buffer_[isp].GetRealData(DiagIdx::w).size());
-        const int total_reducedDiag_size = copy_size + init_size;
-
-        // allocate array size for reduced diagnostic buffer array
-        m_particles_buffer_[isp].resize(total_reducedDiag_size);
-
-        // Data pointers to particle attributes //
-        ParticleReal* const AMREX_RESTRICT wp_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::w).data();
-        ParticleReal* const AMREX_RESTRICT x_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::x).data();
-        ParticleReal* const AMREX_RESTRICT y_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::y).data();
-        ParticleReal* const AMREX_RESTRICT z_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::z).data();
-        ParticleReal* const AMREX_RESTRICT ux_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::ux).data();
-        ParticleReal* const AMREX_RESTRICT uy_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::uy).data();
-        ParticleReal* const AMREX_RESTRICT uz_buff =
-              m_particles_buffer_[isp].GetRealData(DiagIdx::uz).data();
-
-        // Selective copy of particle data from tmp array to reduced buffer
-        // array on the GPU using the flag value and index location.
-        amrex::ParallelFor(np,
-        [=] AMREX_GPU_DEVICE(int i)
-        {
-            if (Flag[i] == 1)
-            {
-               const int loc = IndexLocation[i] + init_size;
-               wp_buff[loc] = wp_temp[i];
-               x_buff[loc]  = x_temp[i];
-               y_buff[loc]  = y_temp[i];
-               z_buff[loc]  = z_temp[i];
-               ux_buff[loc] = ux_temp[i];
-               uy_buff[loc] = uy_temp[i];
-               uz_buff[loc] = uz_temp[i];
-            }
-        });
-
-    }
-}
diff --git a/Source/Diagnostics/BackTransformedDiagnostic_fwd.H b/Source/Diagnostics/BackTransformedDiagnostic_fwd.H
deleted file mode 100644
index 2766d2e27..000000000
--- a/Source/Diagnostics/BackTransformedDiagnostic_fwd.H
+++ /dev/null
@@ -1,15 +0,0 @@
-/* Copyright 2021 Luca Fedeli, Axel Huebl
- *
- * This file is part of WarpX.
- *
- * License: BSD-3-Clause-LBNL
- */
-
-#ifndef WARPX_BACK_TRANSFORMED_DIAGNOSTICS_FWD_H
-#define WARPX_BACK_TRANSFORMED_DIAGNOSTICS_FWD_H
-
-class LabFrameDiag;
-class LabFrameSnapShot;
-class BackTransformedDiagnostic;
-
-#endif /* WARPX_BACK_TRANSFORMED_DIAGNOSTICS_FWD_H */
diff --git a/Source/Diagnostics/CMakeLists.txt b/Source/Diagnostics/CMakeLists.txt
index 4df0f68f8..533583394 100644
--- a/Source/Diagnostics/CMakeLists.txt
+++ b/Source/Diagnostics/CMakeLists.txt
@@ -1,6 +1,5 @@
 target_sources(WarpX
   PRIVATE
-    BackTransformedDiagnostic.cpp
     Diagnostics.cpp
     FieldIO.cpp
     FullDiagnostics.cpp
diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package
index f8bf69f98..a09b1eaed 100644
--- a/Source/Diagnostics/Make.package
+++ b/Source/Diagnostics/Make.package
@@ -2,7 +2,6 @@ CEXE_sources += MultiDiagnostics.cpp
 CEXE_sources += Diagnostics.cpp
 CEXE_sources += FullDiagnostics.cpp
 CEXE_sources += WarpXIO.cpp
-CEXE_sources += BackTransformedDiagnostic.cpp
 CEXE_sources += ParticleIO.cpp
 CEXE_sources += FieldIO.cpp
 CEXE_sources += SliceDiagnostic.cpp
diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp
index 548c95cd8..1785a04c1 100644
--- a/Source/Diagnostics/WarpXIO.cpp
+++ b/Source/Diagnostics/WarpXIO.cpp
@@ -352,43 +352,3 @@ WarpX::InitFromCheckpoint ()
     mypc->Restart(restart_chkfile);
 
 }
-
-
-std::unique_ptr<MultiFab>
-WarpX::GetCellCenteredData() {
-
-    WARPX_PROFILE("WarpX::GetCellCenteredData()");
-
-    const amrex::IntVect ng(1);
-    const int nc = 10;
-
-    Vector<std::unique_ptr<MultiFab> > cc(finest_level+1);
-
-    for (int lev = 0; lev <= finest_level; ++lev)
-    {
-        cc[lev] = std::make_unique<MultiFab>(grids[lev], dmap[lev], nc, ng );
-
-        int dcomp = 0;
-        // first the electric field
-        AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dmap[lev], dcomp, ng );
-        dcomp += 3;
-        // then the magnetic field
-        AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dmap[lev], dcomp, ng );
-        dcomp += 3;
-        // then the current density
-        AverageAndPackVectorField( *cc[lev], current_fp[lev], dmap[lev], dcomp, ng );
-        dcomp += 3;
-        // then the charge density
-        const std::unique_ptr<MultiFab>& charge_density = mypc->GetChargeDensity(lev);
-        AverageAndPackScalarField( *cc[lev], *charge_density, dmap[lev], dcomp, ng );
-
-        ablastr::utils::communication::FillBoundary(*cc[lev], WarpX::do_single_precision_comms, geom[lev].periodicity());
-    }
-
-    for (int lev = finest_level; lev > 0; --lev)
-    {
-        CoarsenIO::Coarsen( *cc[lev-1], *cc[lev], 0, 0, nc, 0, refRatio(lev-1) );
-    }
-
-    return std::move(cc[0]);
-}
diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp
index 6e04adc3b..42042b7c7 100644
--- a/Source/Evolve/WarpXEvolve.cpp
+++ b/Source/Evolve/WarpXEvolve.cpp
@@ -11,7 +11,6 @@
 #include "WarpX.H"
 
 #include "BoundaryConditions/PML.H"
-#include "Diagnostics/BackTransformedDiagnostic.H"
 #include "Diagnostics/MultiDiagnostics.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
 #include "Evolve/WarpXDtType.H"
@@ -247,15 +246,6 @@ WarpX::Evolve (int numsteps)
 
         ShiftGalileanBoundary();
 
-        if (do_back_transformed_diagnostics) {
-            std::unique_ptr<MultiFab> cell_centered_data = nullptr;
-            if (WarpX::do_back_transformed_fields) {
-                cell_centered_data = GetCellCenteredData();
-            }
-            myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]);
-        }
-
-
         // sync up time
         for (int i = 0; i <= max_level; ++i) {
             t_new[i] = cur_time;
@@ -370,10 +360,6 @@ WarpX::Evolve (int numsteps)
         // End loop on time steps
     }
     multi_diags->FilterComputePackFlushLastTimestep( istep[0] );
-
-    if (do_back_transformed_diagnostics) {
-        myBFD->Flush(geom[0]);
-    }
 }
 
 /* /brief Perform one PIC iteration, without subcycling
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 987bab9d4..3a5e9ed24 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -14,7 +14,6 @@
 #if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
 #   include "BoundaryConditions/PML_RZ.H"
 #endif
-#include "Diagnostics/BackTransformedDiagnostic.H"
 #include "Diagnostics/MultiDiagnostics.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
@@ -438,28 +437,6 @@ WarpX::InitData ()
 void
 WarpX::InitDiagnostics () {
     multi_diags->InitData();
-    if (do_back_transformed_diagnostics) {
-        const Real* current_lo = geom[0].ProbLo();
-        const Real* current_hi = geom[0].ProbHi();
-        Real dt_boost = dt[0];
-        Real boosted_moving_window_v = (moving_window_v - beta_boost*PhysConst::c)/(1 - beta_boost*moving_window_v/PhysConst::c);
-        // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0
-        Real zmin_lab = static_cast<Real>(
-            (current_lo[moving_window_dir] - boosted_moving_window_v*t_new[0])/( (1.+beta_boost)*gamma_boost ));
-        Real zmax_lab = static_cast<Real>(
-            (current_hi[moving_window_dir] - boosted_moving_window_v*t_new[0])/( (1.+beta_boost)*gamma_boost ));
-        myBFD = std::make_unique<BackTransformedDiagnostic>(
-                                               zmin_lab,
-                                               zmax_lab,
-                                               moving_window_v, dt_snapshots_lab,
-                                               num_snapshots_lab,
-                                               dt_slice_snapshots_lab,
-                                               num_slice_snapshots_lab,
-                                               gamma_boost, t_new[0], dt_boost,
-                                               moving_window_dir, geom[0],
-                                               slice_realbox,
-                                               particle_slice_width_lab);
-    }
     reduced_diags->InitData();
 }
 
diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp
index b810723fe..001fc4631 100644
--- a/Source/Particles/LaserParticleContainer.cpp
+++ b/Source/Particles/LaserParticleContainer.cpp
@@ -85,7 +85,6 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
 {
     charge = 1.0;
     mass = std::numeric_limits<Real>::max();
-    do_back_transformed_diagnostics = 0;
 
     ParmParse pp_laser_name(m_laser_name);
 
diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H
index 7f355810e..59815d150 100644
--- a/Source/Particles/MultiParticleContainer.H
+++ b/Source/Particles/MultiParticleContainer.H
@@ -236,9 +236,6 @@ public:
 
     int nSpecies() const {return species_names.size();}
 
-    int nSpeciesBackTransformedDiagnostics() const {return nspecies_back_transformed_diagnostics;}
-    int mapSpeciesBackTransformedDiagnostics(int i) const {return map_species_back_transformed_diagnostics[i];}
-    int doBackTransformedDiagnostics() const {return do_back_transformed_diagnostics;}
     /** Whether back-transformed diagnostics need to be performed for any plasma species.
      *
      * \param[in] do_back_transformed_particles The parameter to set if back-transformed particles are set to true/false
@@ -263,12 +260,6 @@ public:
         return std::count( v.begin(), v.end(), fromMainGrid );
     }
 
-    void GetLabFrameData(const std::string& snapshot_name,
-                         const int i_lab, const int direction,
-                         const amrex::Real z_old, const amrex::Real z_new,
-                         const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt,
-                         amrex::Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const;
-
     // Inject particles during the simulation (for particles entering the
     // simulation domain after some iterations, due to flowing plasma and/or
     // moving window).
@@ -494,12 +485,6 @@ private:
 
     void mapSpeciesProduct ();
 
-    // Number of species dumped in BackTransformedDiagnostics
-    int nspecies_back_transformed_diagnostics = 0;
-    // map_species_back_transformed_diagnostics[i] is the species ID in
-    // MultiParticleContainer for 0<i<nspecies_back_transformed_diagnostics
-    std::vector<int> map_species_back_transformed_diagnostics;
-    int do_back_transformed_diagnostics = 0;
     bool m_do_back_transformed_particles = false;
 
     void MFItInfoCheckTiling(const WarpXParticleContainer& /*pc_src*/) const noexcept
diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 7200df69a..9cbe898cd 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -119,20 +119,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core)
 
     pc_tmp = std::make_unique<PhysicalParticleContainer>(amr_core);
 
-    // Compute the number of species for which lab-frame data is dumped
-    // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer
-    // particle IDs in map_species_lab_diags.
-    map_species_back_transformed_diagnostics.resize(nspecies);
-    nspecies_back_transformed_diagnostics = 0;
-    for (int i=0; i<nspecies; i++){
-        auto& pc = allcontainers[i];
-        if (pc->do_back_transformed_diagnostics){
-            map_species_back_transformed_diagnostics[nspecies_back_transformed_diagnostics] = i;
-            do_back_transformed_diagnostics = 1;
-            nspecies_back_transformed_diagnostics += 1;
-        }
-    }
-
     // Setup particle collisions
     collisionhandler = std::make_unique<CollisionHandler>(this);
 
@@ -726,68 +712,6 @@ MultiParticleContainer::SetParticleDistributionMap (int lev, DistributionMapping
     }
 }
 
-void
-MultiParticleContainer
-::GetLabFrameData (const std::string& /*snapshot_name*/,
-                   const int /*i_lab*/, const int direction,
-                   const Real z_old, const Real z_new,
-                   const Real t_boost, const Real t_lab, const Real dt,
-                   Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const
-{
-
-    WARPX_PROFILE("MultiParticleContainer::GetLabFrameData()");
-
-    // Loop over particle species
-    for (int i = 0; i < nspecies_back_transformed_diagnostics; ++i){
-        int isp = map_species_back_transformed_diagnostics[i];
-        WarpXParticleContainer* pc = allcontainers[isp].get();
-        WarpXParticleContainer::DiagnosticParticles diagnostic_particles;
-        pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles);
-        // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData
-        // where "lev" is the AMR level and "index" is a [grid index][tile index] pair.
-
-        // Loop over AMR levels
-        for (int lev = 0; lev <= pc->finestLevel(); ++lev){
-            // Loop over [grid index][tile index] pairs
-            // and Fills parts[species number i] with particle data from all grids and
-            // tiles in diagnostic_particles. parts contains particles from all
-            // AMR levels indistinctly.
-            for (const auto& dp : diagnostic_particles[lev]){
-                // it->first is the [grid index][tile index] key
-                // it->second is the corresponding
-                // WarpXParticleContainer::DiagnosticParticleData value
-                parts[i].GetRealData(DiagIdx::w).insert(  parts[i].GetRealData(DiagIdx::w  ).end(),
-                                                          dp.second.GetRealData(DiagIdx::w  ).begin(),
-                                                          dp.second.GetRealData(DiagIdx::w  ).end());
-
-                parts[i].GetRealData(DiagIdx::x).insert(  parts[i].GetRealData(DiagIdx::x  ).end(),
-                                                          dp.second.GetRealData(DiagIdx::x  ).begin(),
-                                                          dp.second.GetRealData(DiagIdx::x  ).end());
-
-                parts[i].GetRealData(DiagIdx::y).insert(  parts[i].GetRealData(DiagIdx::y  ).end(),
-                                                          dp.second.GetRealData(DiagIdx::y  ).begin(),
-                                                          dp.second.GetRealData(DiagIdx::y  ).end());
-
-                parts[i].GetRealData(DiagIdx::z).insert(  parts[i].GetRealData(DiagIdx::z  ).end(),
-                                                          dp.second.GetRealData(DiagIdx::z  ).begin(),
-                                                          dp.second.GetRealData(DiagIdx::z  ).end());
-
-                parts[i].GetRealData(DiagIdx::ux).insert(  parts[i].GetRealData(DiagIdx::ux).end(),
-                                                           dp.second.GetRealData(DiagIdx::ux).begin(),
-                                                           dp.second.GetRealData(DiagIdx::ux).end());
-
-                parts[i].GetRealData(DiagIdx::uy).insert(  parts[i].GetRealData(DiagIdx::uy).end(),
-                                                           dp.second.GetRealData(DiagIdx::uy).begin(),
-                                                           dp.second.GetRealData(DiagIdx::uy).end());
-
-                parts[i].GetRealData(DiagIdx::uz).insert(  parts[i].GetRealData(DiagIdx::uz).end(),
-                                                           dp.second.GetRealData(DiagIdx::uz).begin(),
-                                                           dp.second.GetRealData(DiagIdx::uz).end());
-            }
-        }
-    }
-}
-
 /* \brief Continuous injection for particles initially outside of the domain.
  * \param injection_box: Domain where new particles should be injected.
  * Loop over all WarpXParticleContainer in MultiParticleContainer and
diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp
index 29a835735..2de928aed 100644
--- a/Source/Particles/PhotonParticleContainer.cpp
+++ b/Source/Particles/PhotonParticleContainer.cpp
@@ -127,8 +127,7 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti,
 #endif
 
     auto copyAttribs = CopyParticleAttribs(pti, tmp_particle_data, offset);
-    int do_copy = (WarpX::do_back_transformed_diagnostics &&
-                   do_back_transformed_diagnostics && a_dt_type!=DtType::SecondHalf);
+    int do_copy = (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf) );
 
     const auto GetPosition = GetParticlePosition(pti, offset);
     auto SetPosition = SetParticlePosition(pti, offset);
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index 163658102..2f52688fb 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -240,12 +240,6 @@ public:
                     const int n_external_attr_int,
                     const amrex::RandomEngine& engine) override final;
 
-    virtual void GetParticleSlice (
-        const int direction, const amrex::Real z_old,
-        const amrex::Real z_new, const amrex::Real t_boost,
-        const amrex::Real t_lab, const amrex::Real dt,
-        DiagnosticParticles& diagnostic_particles) final;
-
 /**
  * \brief Apply NCI Godfrey filter to all components of E and B before gather
  * \param lev MR level
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index e79d526f1..b5ea417c0 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -264,9 +264,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
     utils::parser::queryWithParser(
         pp_species_name, "self_fields_max_iters", self_fields_max_iters);
     pp_species_name.query("self_fields_verbosity", self_fields_verbosity);
-    // Whether to plot back-transformed (lab-frame) diagnostics
-    // for this species.
-    pp_species_name.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics);
 
     pp_species_name.query("do_field_ionization", do_field_ionization);
 
@@ -1862,8 +1859,7 @@ PhysicalParticleContainer::Evolve (int lev,
 
     bool has_buffer = cEx || cjx;
 
-    if ( (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) ||
-         (m_do_back_transformed_particles) )
+    if (m_do_back_transformed_particles)
     {
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
@@ -2478,191 +2474,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
     }
 }
 
-void
-PhysicalParticleContainer::GetParticleSlice (
-    const int direction, const Real z_old,
-    const Real z_new, const Real t_boost,
-    const Real t_lab, const Real dt,
-    DiagnosticParticles& diagnostic_particles)
-{
-    WARPX_PROFILE("PhysicalParticleContainer::GetParticleSlice()");
-
-    // Assume that the boost in the positive z direction.
-#if defined(WARPX_DIM_1D_Z)
-    AMREX_ALWAYS_ASSERT(direction == 0);
-#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-    AMREX_ALWAYS_ASSERT(direction == 1);
-#else
-    AMREX_ALWAYS_ASSERT(direction == 2);
-#endif
-
-    // Note the the slice should always move in the negative boost direction.
-    AMREX_ALWAYS_ASSERT(z_new < z_old);
-
-    AMREX_ALWAYS_ASSERT(do_back_transformed_diagnostics == 1);
-
-    const int nlevs = std::max(0, finestLevel()+1);
-    diagnostic_particles.resize(finestLevel()+1);
-
-    for (int lev = 0; lev < nlevs; ++lev) {
-        // first we touch each map entry in serial
-        for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
-        {
-            auto index = std::make_pair(pti.index(), pti.LocalTileIndex());
-            diagnostic_particles[lev][index];
-        }
-
-#ifdef AMREX_USE_OMP
-#pragma omp parallel
-#endif
-        {
-            // Temporary arrays to store copy_flag and copy_index
-            // for particles that cross the z-slice
-            // These arrays are defined before the WarpXParIter to prevent them
-            // from going out of scope after each iteration, while the kernels
-            // may still need access to them.
-            // Note that the destructor for WarpXParIter is synchronized.
-            amrex::Gpu::DeviceVector<int> FlagForPartCopy;
-            amrex::Gpu::DeviceVector<int> IndexForPartCopy;
-            for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
-            {
-                auto index = std::make_pair(pti.index(), pti.LocalTileIndex());
-
-                const auto GetPosition = GetParticlePosition(pti);
-
-                auto& attribs = pti.GetAttribs();
-                ParticleReal* const AMREX_RESTRICT wpnew = attribs[PIdx::w].dataPtr();
-                ParticleReal* const AMREX_RESTRICT uxpnew = attribs[PIdx::ux].dataPtr();
-                ParticleReal* const AMREX_RESTRICT uypnew = attribs[PIdx::uy].dataPtr();
-                ParticleReal* const AMREX_RESTRICT uzpnew = attribs[PIdx::uz].dataPtr();
-
-                ParticleReal* const AMREX_RESTRICT
-                  xpold = tmp_particle_data[lev][index][TmpIdx::xold].dataPtr();
-                ParticleReal* const AMREX_RESTRICT
-                  ypold = tmp_particle_data[lev][index][TmpIdx::yold].dataPtr();
-                ParticleReal* const AMREX_RESTRICT
-                  zpold = tmp_particle_data[lev][index][TmpIdx::zold].dataPtr();
-                ParticleReal* const AMREX_RESTRICT
-                  uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr();
-                ParticleReal* const AMREX_RESTRICT
-                  uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr();
-                ParticleReal* const AMREX_RESTRICT
-                  uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr();
-
-                const long np = pti.numParticles();
-
-                Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c;
-                Real inv_c2 = 1.0_rt/PhysConst::c/PhysConst::c;
-
-                FlagForPartCopy.resize(np);
-                IndexForPartCopy.resize(np);
-
-                int* const AMREX_RESTRICT Flag = FlagForPartCopy.dataPtr();
-                int* const AMREX_RESTRICT IndexLocation = IndexForPartCopy.dataPtr();
-
-                //Flag particles that need to be copied if they cross the z_slice
-                amrex::ParallelFor(np,
-                [=] AMREX_GPU_DEVICE(int i)
-                {
-                    ParticleReal xp, yp, zp;
-                    GetPosition(i, xp, yp, zp);
-                    Flag[i] = 0;
-                    if ( (((zp >= z_new) && (zpold[i] <= z_old)) ||
-                          ((zp <= z_new) && (zpold[i] >= z_old))) )
-                    {
-                        Flag[i] = 1;
-                    }
-                });
-
-                // exclusive scan to obtain location indices using flag values
-                // These location indices are used to copy data from
-                // src to dst when the copy-flag is set to 1.
-                const int total_partdiag_size = amrex::Scan::ExclusiveSum(np,Flag,IndexLocation);
-
-                // allocate array size for diagnostic particle array
-                diagnostic_particles[lev][index].resize(total_partdiag_size);
-
-                amrex::Real gammaboost = WarpX::gamma_boost;
-                amrex::Real betaboost = WarpX::beta_boost;
-                amrex::Real Phys_c = PhysConst::c;
-
-                ParticleReal* const AMREX_RESTRICT diag_wp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::w).data();
-                ParticleReal* const AMREX_RESTRICT diag_xp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::x).data();
-                ParticleReal* const AMREX_RESTRICT diag_yp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::y).data();
-                ParticleReal* const AMREX_RESTRICT diag_zp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::z).data();
-                ParticleReal* const AMREX_RESTRICT diag_uxp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::ux).data();
-                ParticleReal* const AMREX_RESTRICT diag_uyp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::uy).data();
-                ParticleReal* const AMREX_RESTRICT diag_uzp =
-                diagnostic_particles[lev][index].GetRealData(DiagIdx::uz).data();
-
-                // Copy particle data to diagnostic particle array on the GPU
-                //  using flag and index values
-                amrex::ParallelFor(np,
-                [=] AMREX_GPU_DEVICE(int i)
-                {
-                    ParticleReal xp_new, yp_new, zp_new;
-                    GetPosition(i, xp_new, yp_new, zp_new);
-                    if (Flag[i] == 1)
-                    {
-                         // Lorentz Transform particles to lab-frame
-                         const Real gamma_new_p = std::sqrt(1.0_rt + inv_c2*
-                                                  (uxpnew[i]*uxpnew[i]
-                                                 + uypnew[i]*uypnew[i]
-                                                 + uzpnew[i]*uzpnew[i]));
-                         const Real t_new_p = gammaboost*t_boost - uzfrm*zp_new*inv_c2;
-                         const Real z_new_p = gammaboost*(zp_new + betaboost*Phys_c*t_boost);
-                         const Real uz_new_p = gammaboost*uzpnew[i] - gamma_new_p*uzfrm;
-
-                         const Real gamma_old_p = std::sqrt(1.0_rt + inv_c2*
-                                                  (uxpold[i]*uxpold[i]
-                                                 + uypold[i]*uypold[i]
-                                                 + uzpold[i]*uzpold[i]));
-                         const Real t_old_p = gammaboost*(t_boost - dt)
-                                              - uzfrm*zpold[i]*inv_c2;
-                         const Real z_old_p = gammaboost*(zpold[i]
-                                              + betaboost*Phys_c*(t_boost-dt));
-                         const Real uz_old_p = gammaboost*uzpold[i]
-                                              - gamma_old_p*uzfrm;
-
-                         // interpolate in time to t_lab
-                         const Real weight_old = (t_new_p - t_lab)
-                                               / (t_new_p - t_old_p);
-                         const Real weight_new = (t_lab - t_old_p)
-                                               / (t_new_p - t_old_p);
-
-                         const Real xp = xpold[i]*weight_old + xp_new*weight_new;
-                         const Real yp = ypold[i]*weight_old + yp_new*weight_new;
-                         const Real zp = z_old_p*weight_old  + z_new_p*weight_new;
-
-                         const Real uxp = uxpold[i]*weight_old
-                                        + uxpnew[i]*weight_new;
-                         const Real uyp = uypold[i]*weight_old
-                                        + uypnew[i]*weight_new;
-                         const Real uzp = uz_old_p*weight_old
-                                        + uz_new_p  *weight_new;
-
-                         const int loc = IndexLocation[i];
-                         diag_wp[loc] = wpnew[i];
-                         diag_xp[loc] = xp;
-                         diag_yp[loc] = yp;
-                         diag_zp[loc] = zp;
-                         diag_uxp[loc] = uxp;
-                         diag_uyp[loc] = uyp;
-                         diag_uzp[loc] = uzp;
-                    }
-                });
-                Gpu::synchronize(); // because of FlagForPartCopy & IndexForPartCopy
-            }
-        }
-    }
-}
-
 /* \brief Inject particles during the simulation
  * \param injection_box: domain where particles should be injected.
  */
@@ -2766,10 +2577,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti,
     ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr() + offset;
     ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr() + offset;
 
-    int do_copy = ( (WarpX::do_back_transformed_diagnostics
-                     && do_back_transformed_diagnostics
-                     && a_dt_type!=DtType::SecondHalf)
-                  || (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf)) );
+    int do_copy = (m_do_back_transformed_particles && (a_dt_type!=DtType::SecondHalf) );
     CopyParticleAttribs copyAttribs;
     if (do_copy) {
         copyAttribs = CopyParticleAttribs(pti, tmp_particle_data, offset);
diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H
index 6f0aaf8e8..cb5542025 100644
--- a/Source/Particles/WarpXParticleContainer.H
+++ b/Source/Particles/WarpXParticleContainer.H
@@ -137,11 +137,6 @@ public:
 
     virtual void PostRestart () = 0;
 
-    virtual void GetParticleSlice(const int /*direction*/, const amrex::Real /*z_old*/,
-                                  const amrex::Real /*z_new*/, const amrex::Real /*t_boost*/,
-                                  const amrex::Real /*t_lab*/, const amrex::Real /*dt*/,
-                                  DiagnosticParticles& /*diagnostic_particles*/) {}
-
     void AllocData ();
 
     /**
@@ -315,7 +310,6 @@ public:
     // split along diagonals (0) or axes (1)
     int split_type = 0;
 
-    int doBackTransformedDiagnostics () const { return do_back_transformed_diagnostics; }
     /** Whether back-transformed diagnostics need to be performed for a particular species.
      *
      * \param[in] do_back_transformed_particles The parameter to set if back-transformed particles are set to true/false
@@ -408,7 +402,6 @@ protected:
 
     int do_resampling = 0;
 
-    int do_back_transformed_diagnostics = 1;
     /** Whether back-transformed diagnostics is turned on for the corresponding species.*/
     bool m_do_back_transformed_particles = false;
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 1cad22d40..ecedba62b 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -13,7 +13,6 @@
 #define WARPX_H_
 
 #include "BoundaryConditions/PML_fwd.H"
-#include "Diagnostics/BackTransformedDiagnostic_fwd.H"
 #include "Diagnostics/MultiDiagnostics_fwd.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags_fwd.H"
 #include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H"
@@ -273,22 +272,6 @@ public:
     //! If true, the initial conditions from random number generators are serialized (useful for reproducible testing with OpenMP)
     static bool serialize_initial_conditions;
 
-    //! If true, then lab-frame data will be computed for boosted frame simulations
-    //! with customized output format
-    static bool do_back_transformed_diagnostics;
-    //! Name of the back-transformed diagnostics output directory
-    static std::string lab_data_directory;
-    //! Number of back-tranformed snapshots in the lab-frame
-    static int  num_snapshots_lab;
-    //! Time interval in lab-frame between the back-transformed snapshots
-    static amrex::Real dt_snapshots_lab;
-    //! If true, then lab-frame data will be computed for the fields and flushed out
-    //! in customized format
-    static bool do_back_transformed_fields;
-    //! If true, then lab-frame data will be computed for the particles and flushed out
-    //! in customized format
-    static bool do_back_transformed_particles;
-
     //! Lorentz factor of the boosted frame in which a boosted-frame simulation is run
     static amrex::Real gamma_boost;
     //! Beta value corresponding to the Lorentz factor of the boosted frame of the simulation
@@ -758,12 +741,6 @@ public:
     static amrex::Real moving_window_v;
     static bool fft_do_time_averaging;
 
-    // slice generation //
-    static int num_slice_snapshots_lab;
-    static amrex::Real dt_slice_snapshots_lab;
-    static amrex::Real particle_slice_width_lab;
-    amrex::RealBox getSliceRealBox() const {return slice_realbox;}
-
     // these should be private, but can't due to Cuda limitations
     static void ComputeDivB (amrex::MultiFab& divB, int const dcomp,
                              const std::array<const amrex::MultiFab* const, 3>& B,
@@ -1077,8 +1054,6 @@ private:
     /** Check the requested resources and write performance hints */
     void PerformanceHints ();
 
-    std::unique_ptr<amrex::MultiFab> GetCellCenteredData();
-
     void BuildBufferMasks ();
     void BuildBufferMasksInBox ( const amrex::Box tbx, amrex::IArrayBox &buffer_mask,
                                  const amrex::IArrayBox &guard_mask, const int ng );
@@ -1151,9 +1126,6 @@ private:
     std::unique_ptr<MultiParticleContainer> mypc;
     std::unique_ptr<MultiDiagnostics> multi_diags;
 
-    // Boosted Frame Diagnostics
-    std::unique_ptr<BackTransformedDiagnostic> myBFD;
-
     //
     // Fields: First array for level, second for direction
     //
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 2e0fe7c3e..9ee07811d 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -12,7 +12,6 @@
 #include "WarpX.H"
 
 #include "BoundaryConditions/PML.H"
-#include "Diagnostics/BackTransformedDiagnostic.H"
 #include "Diagnostics/MultiDiagnostics.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags.H"
 #include "EmbeddedBoundary/WarpXFaceInfoBox.H"
@@ -169,17 +168,6 @@ int WarpX::num_mirrors = 0;
 utils::parser::IntervalsParser WarpX::sort_intervals;
 amrex::IntVect WarpX::sort_bin_size(AMREX_D_DECL(1,1,1));
 
-bool WarpX::do_back_transformed_diagnostics = false;
-std::string WarpX::lab_data_directory = "lab_frame_data";
-int  WarpX::num_snapshots_lab = std::numeric_limits<int>::lowest();
-Real WarpX::dt_snapshots_lab  = std::numeric_limits<Real>::lowest();
-bool WarpX::do_back_transformed_fields = true;
-bool WarpX::do_back_transformed_particles = true;
-
-int  WarpX::num_slice_snapshots_lab = 0;
-Real WarpX::dt_slice_snapshots_lab;
-Real WarpX::particle_slice_width_lab = 0.0_rt;
-
 bool WarpX::do_dynamic_scheduling = true;
 
 int WarpX::electrostatic_solver_id;
@@ -270,7 +258,6 @@ WarpX::WarpX ()
             current_injection_position = geom[0].ProbLo(moving_window_dir);
         }
     }
-    do_back_transformed_particles = mypc->doBackTransformedDiagnostics();
 
     // Particle Boundary Buffer (i.e., scraped particles on boundary)
     m_particle_boundary_buffer = std::make_unique<ParticleBoundaryBuffer>();
@@ -644,49 +631,6 @@ WarpX::ReadParameters ()
             moving_window_v *= PhysConst::c;
         }
 
-        pp_warpx.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics);
-        if (do_back_transformed_diagnostics) {
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0,
-                   "gamma_boost must be > 1 to use the boosted frame diagnostic.");
-
-            pp_warpx.query("lab_data_directory", lab_data_directory);
-
-            std::string s;
-            pp_warpx.get("boost_direction", s);
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
-                   "The boosted frame diagnostic currently only works if the boost is in the z direction.");
-
-            utils::parser::queryWithParser(
-                pp_warpx, "num_snapshots_lab", num_snapshots_lab);
-
-            // Read either dz_snapshots_lab or dt_snapshots_lab
-            Real dz_snapshots_lab = 0;
-            const bool dt_snapshots_specified =
-                utils::parser::queryWithParser(pp_warpx, "dt_snapshots_lab", dt_snapshots_lab);
-            const bool dz_snapshots_specified =
-                utils::parser::queryWithParser(pp_warpx, "dz_snapshots_lab", dz_snapshots_lab);
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-                dt_snapshots_specified || dz_snapshots_specified,
-                "When using back-transformed diagnostics, user should specify either dz_snapshots_lab or dt_snapshots_lab.");
-
-            if (dz_snapshots_specified){
-                dt_snapshots_lab = dz_snapshots_lab/PhysConst::c;
-            }
-
-            utils::parser::getWithParser(pp_warpx, "gamma_boost", gamma_boost);
-
-            pp_warpx.query("do_back_transformed_fields", do_back_transformed_fields);
-
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window,
-                   "The moving window should be on if using the boosted frame diagnostic.");
-
-            pp_warpx.get("moving_window_dir", s);
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"),
-                   "The boosted frame diagnostic currently only works if the moving window is in the z direction.");
-        }
-
         electrostatic_solver_id = GetAlgorithmInteger(pp_warpx, "do_electrostatic");
         // if an electrostatic solver is used, set the Maxwell solver to None
         if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) {
@@ -1451,19 +1395,6 @@ WarpX::ReadParameters ()
                 slice_cr_ratio[idim] = slice_crse_ratio[idim];
             }
         }
-
-        if (do_back_transformed_diagnostics) {
-            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0,
-                "gamma_boost must be > 1 to use the boost frame diagnostic");
-            utils::parser::queryWithParser(
-                pp_slice, "num_slice_snapshots_lab", num_slice_snapshots_lab);
-            if (num_slice_snapshots_lab > 0) {
-                utils::parser::getWithParser(
-                    pp_slice, "dt_slice_snapshots_lab", dt_slice_snapshots_lab );
-                utils::parser::getWithParser(
-                    pp_slice, "particle_slice_width_lab",particle_slice_width_lab);
-            }
-        }
     }
 }
 
@@ -1571,6 +1502,68 @@ WarpX::BackwardCompatibility ()
         "Please use the renamed option warpx.serialize_initial_conditions instead."
     );
 
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("do_back_transformed_diagnostics", backward_int),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("lab_data_directory", backward_str),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("num_snapshots_lab", backward_int),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("dt_snapshots_lab", backward_Real),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("dz_snapshots_lab", backward_Real),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("do_back_transformed_fields", backward_int),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_warpx.query("buffer_size", backward_int),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    ParmParse pp_slice("slice");
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_slice.query("num_slice_snapshots_lab", backward_int),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_slice.query("dt_slice_snapshots_lab", backward_Real),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_slice.query("particle_slice_width_lab", backward_Real),
+        "Legacy back-transformed diagnostics are not supported anymore. "
+        "Please use the new syntax for back-transformed diagnostics, see documentation."
+    );
+
     ParmParse pp_interpolation("interpolation");
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
         !pp_interpolation.query("nox", backward_int) &&
diff --git a/Tools/LibEnsemble/sim/inputs b/Tools/LibEnsemble/sim/inputs
index 11e9bb1e7..03e9968e5 100644
--- a/Tools/LibEnsemble/sim/inputs
+++ b/Tools/LibEnsemble/sim/inputs
@@ -48,10 +48,6 @@ warpx.moving_window_v = 1.0 # in units of the speed of light
 #################################
 warpx.gamma_boost = 30.0
 warpx.boost_direction = z
-warpx.do_back_transformed_diagnostics=0 # 1
-warpx.do_back_transformed_fields=0
-warpx.num_snapshots_lab=20
-warpx.dt_snapshots_lab=1.6678204759907604e-10
 
 #################################
 ############ PLASMA #############
@@ -77,7 +73,6 @@ electrons.predefined_profile_name = "parabolic_channel"
 #         predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
 electrons.predefined_profile_params = 0.0       .02       .297      .003        40.e-6   1.7e23
 electrons.do_continuous_injection = 1
-electrons.do_back_transformed_diagnostics=0
 
 ions.charge = q_e
 ions.mass = m_p
@@ -95,7 +90,6 @@ ions.predefined_profile_name = "parabolic_channel"
 #    predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
 ions.predefined_profile_params = 0.0       .02       .297      .003        40.e-6   1.7e23
 ions.do_continuous_injection = 1
-ions.do_back_transformed_diagnostics=0
 
 electrons2.charge = -q_e
 electrons2.mass = m_e
@@ -113,7 +107,6 @@ electrons2.predefined_profile_name = "parabolic_channel"
 #         predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
 electrons2.predefined_profile_params = 0.3485       .02       .297       .003        40.e-6   1.7e23
 electrons2.do_continuous_injection = 1
-electrons2.do_back_transformed_diagnostics=0
 
 ions2.charge = q_e
 ions2.mass = m_p
@@ -131,7 +124,6 @@ ions2.predefined_profile_name = "parabolic_channel"
 #    predefined_profile_params = z_start   ramp_up   plateau   ramp_down   rc       n0
 ions2.predefined_profile_params = 0.3485       .02       .297       .003        40.e-6   1.7e23
 ions2.do_continuous_injection = 1
-ions2.do_back_transformed_diagnostics=0
 
 beam.charge = -q_e
 beam.mass = m_e
@@ -156,7 +148,6 @@ beam.uy_th = 0.
 beam.uz_th = 39.138943248532286
 beam.zinject_plane = 0.02
 beam.rigid_advance = true
-beam.do_back_transformed_diagnostics=1
 
 #################################
 ######### Lens Mirror ###########
diff --git a/Tools/PerformanceTests/automated_test_4_labdiags_2ppc b/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
index ecfd42da1..f49d92acf 100644
--- a/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
+++ b/Tools/PerformanceTests/automated_test_4_labdiags_2ppc
@@ -34,11 +34,6 @@ warpx.moving_window_v = 1.0 # in units of the speed of light
 warpx.gamma_boost = 15.
 warpx.boost_direction = z
 
-# Diagnostics
-warpx.do_back_transformed_diagnostics = 1
-warpx.num_snapshots_lab = 20
-warpx.dt_snapshots_lab = 7.0e-14
-
 # Species
 particles.species_names = electrons ions