Coupling WarpX with an ✨improved✨ version of the QED library (#1198)

* Initial work to couple improved QED module to WarpX

* WIP to couple with WarpX the new QED library

* Continuing work to couple the new version of the QED library with WarpX

* progress towards completing coupling with new version of QED library

* WarpX coupled with new version of QED library

* default behavior is to display table generation progress

* some host device functions are now device only

* fixed bug

* bugfixing

* updating tests

* updated test

* updated test

* added initial version of tests (not working)

* added check and updated a comment

* fixed bug

* added inputfiles and analysis script for new BW tests

* test for BW process are ready

* modified test

* make lgtm happy

* removed TABs

* initial work to add QS tests (not working)

* removed old tests

* fixed bug in script

* changed position of evolution of optical depth

* progress with QSR tests

* improved test

* very low energy photons are always eliminated

* added tests to regression suite

* improved test

* improved tests

* removed redundant parameter

* removed trailing white space

* updated documentation

* fix lgtm warnings

* fixed missing check on chi parameter

* fixed missing check on chi parameter & bugfixing

* improved comments

* increased tolerance in tests

* updated units in test

* now test succeds if the error is extremely small

* updated checksums

* fixed bug

* fixed some unused or uninitialized variables warnings

* now using ignore_unused instead of commenting out some variables

* fixed warnings

* partial fix of a test

* fixed test

* fixed test

* added checksums

* fixed tests

* fixed benchmark for qed_schwinger2

* removed checksums for tests which do no exist anymore

* fixed checksums for several qed tests

* fixed checksums for several qed tests

* fixed checksums

* removed unwanted checksum

* fixed checksum

* removed files which should have been deleted

* add some const

* [skip ci] added some docstrings and some const

* Update Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* added some docstrings and some const

* replaced ManagedVectors with DeviceVectors

* Update Source/Particles/ElementaryProcess/QEDInternals/QedWrapperCommons.H

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* added some const

* removed unwanted assert

* updated comment

* changed position of GPU synchronization directive

* Update Docs/source/running_cpp/parameters.rst

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Examples/Modules/qed/quantum_synchrotron/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Examples/Modules/qed/quantum_synchrotron/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Examples/Modules/qed/breit_wheeler/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Examples/Modules/qed/breit_wheeler/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* add do_plot option to some analysis scripts

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* uncomment a line

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* simplified input scripts for BW tests

* simplified input scripts for QS tests

* removed unwanted files

* simplified analysis script

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* reverted modification to schwinger analysis script

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* remove outdated comment

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Source/Particles/MultiParticleContainer.cpp

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* fix warnings

* made test more robust

* reset benchmark for qed_breit_wheeler_2d

* fixed bug in test

* make test more robust

* made test more robust

* Update Examples/Modules/qed/quantum_synchrotron/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update Examples/Modules/qed/quantum_synchrotron/analysis.py

Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>

* Update run_test.sh

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

Co-authored-by: Tools <warpx@lbl.gov>
Co-authored-by: NeilZaim <49716072+NeilZaim@users.noreply.github.com>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

13 files changed, 1579 insertions, 687 deletions

diff --git a/Examples/Modules/qed/breit_wheeler/analysis.py b/Examples/Modules/qed/breit_wheeler/analysis.py
new file mode 100755
index 000000000..280e8e2a8
--- /dev/null
+++ b/Examples/Modules/qed/breit_wheeler/analysis.py
@@ -0,0 +1,298 @@
+#! /usr/bin/env python
+
+# Copyright 2019 Luca Fedeli, Maxence Thevenet
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# -*- coding: utf-8 -*-
+
+import yt
+import numpy as np
+import sys
+import scipy.special as spe
+import scipy.integrate as integ
+import scipy.stats as st
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+import matplotlib.pyplot as plt
+
+# This script performs detailed checks of the Breit-Wheeler pair production process.
+# Four populations of photons are initialized with different momenta in different
+# directions in a background EM field (with non-zero components along each direction).
+# Specifically the script checks that:
+#
+# - The expected number of generated pairs n_pairs is in agreement with theory
+#   (the maximum tolerated error is 5*sqrt(n_pairs)
+# - The weight of the generated particles is equal to the weight of the photon
+# - Momenta of the residual photons are still equal to the original momentum
+# - The generated particles are emitted in the right direction
+# - Total energy is conserved in each event
+# - The energy distribution of the generated particles is in agreement with theory
+# - The optical depths of the product species are correctly initialized (QED effects are
+#   enabled for product species too).
+#
+# More details on the theoretical formulas used in this script can be found in
+# the Jupyter notebook picsar/src/multi_physics/QED_tests/validation/validation.ipynb
+#
+# References:
+# 1) R. Duclous et al 2011 Plasma Phys. Control. Fusion 53 015009
+# 2) A. Gonoskov et al. 2015 Phys. Rev. E 92, 023305
+# 3) M. Lobet. PhD thesis "Effets radiatifs et d'electrodynamique
+#    quantique dans l'interaction laser-matiere ultra-relativiste"
+#    URL: https://tel.archives-ouvertes.fr/tel-01314224
+
+
+# Tolerances
+tol = 1.e-8
+tol_red = 2.e-2
+
+# Physical constants (from CODATA 2018, see: https://physics.nist.gov/cuu/Constants/index.html )
+me = 9.1093837015e-31 #electron mass
+c = 299792458 #speed of light
+hbar = 6.62607015e-34/(2*np.pi) #reduced Plank constant
+fine_structure = 7.2973525693e-3 #fine structure constant
+qe = 1.602176634e-19#elementary charge
+E_s = (me**2 * c**3)/(qe * hbar) #Schwinger E field
+B_s = E_s/c #Schwinger B field
+
+mec = me*c
+mec2 = mec*c
+#______________
+
+# Initial parameters
+spec_names_phot = ["p1", "p2", "p3", "p4"]
+spec_names_ele = ["ele1", "ele2", "ele3", "ele4"]
+spec_names_pos = ["pos1", "pos2", "pos3", "pos4"]
+initial_momenta = [
+    np.array([2000.0,0,0])*mec,
+    np.array([0.0,5000.0,0.0])*mec,
+    np.array([0.0,0.0,10000.0])*mec,
+    np.array([57735.02691896, 57735.02691896, 57735.02691896])*mec
+]
+initial_particle_number = 1048576
+
+E_f = np.array([-2433321316961438., 973328526784575., 1459992790176863.])
+B_f = np.array([2857142.85714286, 4285714.28571428, 8571428.57142857])
+
+NNS = [128,128,128,128] #bins for energy distribution comparison.
+#______________
+
+def calc_chi_gamma(p, E, B):
+    pnorm = np.linalg.norm(p)
+    v = c*(p/pnorm)
+    EpvvecB = E + np.cross(v,B)
+    vdotEoverc = np.dot(v,E)/c
+    ff = np.sqrt(np.dot(EpvvecB,EpvvecB) - np.dot(vdotEoverc,vdotEoverc))
+    gamma_phot = pnorm/mec
+    return gamma_phot*ff/E_s
+
+#Auxiliary functions
+@np.vectorize
+def BW_inner(x):
+    return integ.quad(lambda s: np.sqrt(s)*spe.kv(1./3., 2./3. * s**(3./2.)), x, np.inf)[0]
+
+def BW_X(chi_phot, chi_ele):
+    div = (chi_ele*(chi_phot-chi_ele))
+    div = np.where(np.logical_and(chi_phot > chi_ele, chi_ele != 0), div, 1.0);
+    res = np.where(np.logical_and(chi_phot > chi_ele, chi_ele != 0), np.power(chi_phot/div, 2./3.), np.inf)
+    return res
+
+def BW_F(chi_phot, chi_ele):
+    X = BW_X(chi_phot, chi_ele)
+    res = np.where(np.logical_or(chi_phot == chi_ele, chi_ele == 0), 0.0,
+         BW_inner(X) - (2.0 - chi_phot* X**(3./2.))*spe.kv(2./3., 2./3. * X**(3./2.)) )
+    return res
+
+@np.vectorize
+def BW_T(chi_phot):
+    coeff = 1./(np.pi * np.sqrt(3.) * (chi_phot**2))
+    return coeff*integ.quad(lambda chi_ele: BW_F(chi_phot, chi_ele), 0, chi_phot)[0]
+
+def small_diff(vv, val):
+    if(val != 0.0):
+        return np.max(np.abs((vv - val)/val)) < tol
+    else:
+        return np.max(np.abs(vv)) < tol
+#__________________
+
+# Breit-Wheeler total and differential cross sections
+def BW_dN_dt(chi_phot, gamma_phot):
+    coeff_BW = fine_structure * me*c**2/hbar
+    return coeff_BW*BW_T(chi_phot)*(chi_phot/gamma_phot)
+
+def BW_d2N_dt_dchi(chi_phot, gamma_phot, chi_ele):
+    coeff_BW = fine_structure * me*c**2/hbar
+    return coeff_BW*BW_F(chi_phot, chi_ele)*(gamma_phot/gamma_phot)
+#__________________
+
+# Get data for a species
+def get_spec(ytdata, specname, is_photon):
+    px = ytdata[specname,"particle_momentum_x"].v
+    pz = ytdata[specname,"particle_momentum_z"].v
+    py = ytdata[specname,"particle_momentum_y"].v
+
+    w = ytdata[specname,"particle_weighting"].v
+
+    if (is_photon):
+        opt = ytdata[specname,"particle_optical_depth_BW"].v
+    else:
+        opt = ytdata[specname,"particle_optical_depth_QSR"].v
+
+    return {"px" : px, "py" : py, "pz" : pz, "w" : w, "opt" : opt}
+
+# Individual tests
+
+def check_number_of_pairs(ytdataset, phot_name, ele_name, pos_name, chi_phot, gamma_phot, dt, particle_number):
+    dNBW_dt_theo = BW_dN_dt(chi_phot, gamma_phot)
+    expected_pairs = (1.-np.exp(-dNBW_dt_theo*dt))*particle_number
+    expected_pairs_tolerance = 5.0*np.sqrt(expected_pairs)
+    n_ele = ytdataset.particle_type_counts[ele_name]
+    n_pos = ytdataset.particle_type_counts[pos_name]
+    n_phot = ytdataset.particle_type_counts[phot_name]
+    n_lost = initial_particle_number-n_phot
+    assert((n_ele == n_pos) and (n_ele == n_lost))
+    assert( np.abs(n_ele-expected_pairs) < expected_pairs_tolerance)
+    print("  [OK] generated pair number is within expectations")
+    return n_lost
+
+def check_weights(phot_data, ele_data, pos_data):
+    assert(np.all(phot_data["w"] == phot_data["w"][0]))
+    assert(np.all(ele_data["w"]  == phot_data["w"][0]))
+    assert(np.all(pos_data["w"]  == phot_data["w"][0]))
+    print("  [OK] particles weights are the expected ones")
+
+def check_momenta(phot_data, ele_data, pos_data, p0, p_ele, p_pos):
+    assert(small_diff(phot_data["px"], p0[0]))
+    assert(small_diff(phot_data["py"], p0[1]))
+    assert(small_diff(phot_data["pz"], p0[2]))
+    print("  [OK] residual photons still have initial momentum")
+
+    pdir = p0/np.linalg.norm(p0)
+    assert(small_diff(ele_data["px"]/p_ele, pdir[0]))
+    assert(small_diff(ele_data["py"]/p_ele, pdir[1]))
+    assert(small_diff(ele_data["pz"]/p_ele, pdir[2]))
+    assert(small_diff(pos_data["px"]/p_pos, pdir[0]))
+    assert(small_diff(pos_data["py"]/p_pos, pdir[1]))
+    assert(small_diff(pos_data["pz"]/p_pos, pdir[2]))
+    print("  [OK] pairs move along the initial photon direction")
+
+def check_energy(energy_phot, energy_ele, energy_pos):
+    # Sorting the arrays is required because electrons and positrons are not
+    # necessarily dumped in the same order.
+    s_energy_ele = np.sort(energy_ele)
+    is_energy_pos = np.sort(energy_pos)[::-1]
+    product_energy = s_energy_ele + is_energy_pos
+    assert(small_diff(product_energy, energy_phot))
+    print("  [OK] energy is conserved in each event")
+
+def check_opt_depths(phot_data, ele_data, pos_data):
+    data = (phot_data, ele_data, pos_data)
+    for dd in data:
+        loc, scale = st.expon.fit(dd["opt"])
+        assert( np.abs(loc - 0) < tol_red )
+        assert( np.abs(scale - 1) < tol_red )
+    print("  [OK] optical depth distributions are still exponential")
+
+def check_energy_distrib(energy_ele, energy_pos, gamma_phot,
+        chi_phot, n_lost, NN, idx, do_plot=False):
+    gamma_min = 1.0001
+    gamma_max = gamma_phot-1.0001
+    h_gamma_ele, c_gamma = np.histogram(energy_ele/mec2, bins=NN, range=[gamma_min,gamma_max])
+    h_gamma_pos, _ = np.histogram(energy_pos/mec2, bins=NN, range=[gamma_min,gamma_max])
+
+    cchi_part_min = chi_phot*(gamma_min - 1)/(gamma_phot - 2)
+    cchi_part_max = chi_phot*(gamma_max - 1)/(gamma_phot - 2)
+
+    #Rudimentary integration over npoints for each bin
+    npoints= 20
+    aux_chi = np.linspace(cchi_part_min, cchi_part_max, NN*npoints)
+    distrib = BW_d2N_dt_dchi(chi_phot, gamma_phot, aux_chi)
+    distrib = np.sum(distrib.reshape(-1, npoints),1)
+    distrib = n_lost*distrib/np.sum(distrib)
+
+    if do_plot :
+        # Visual comparison of distributions
+        c_gamma_centered = 0.5*(c_gamma[1:]+c_gamma[:-1])
+        plt.clf()
+        plt.xlabel("γ_particle")
+        plt.ylabel("N")
+        plt.title("χ_photon = {:f}".format(chi_phot))
+        plt.plot(c_gamma_centered, distrib,label="theory")
+        plt.plot(c_gamma_centered, h_gamma_ele,label="BW electrons")
+        plt.plot(c_gamma_centered, h_gamma_pos,label="BW positrons")
+        plt.legend()
+        plt.savefig("case_{:d}".format(idx+1))
+
+    discr_ele = np.abs(h_gamma_ele-distrib)
+    discr_pos = np.abs(h_gamma_pos-distrib)
+    max_discr = 5.0 * np.sqrt(distrib)
+    assert(np.all(discr_ele < max_discr))
+    assert(np.all(discr_pos < max_discr))
+    print("  [OK] energy distribution is within expectations")
+
+#__________________
+
+def check():
+    filename_end = sys.argv[1]
+    data_set_end = yt.load(filename_end)
+
+    sim_time = data_set_end.current_time.to_value()
+    all_data_end = data_set_end.all_data()
+
+    for idx in range(4):
+        phot_name = spec_names_phot[idx]
+        ele_name  = spec_names_ele[idx]
+        pos_name  = spec_names_pos[idx]
+        p0        = initial_momenta[idx]
+
+        p2_phot = p0[0]**2 + p0[1]**2 + p0[2]**2
+        p_phot = np.sqrt(p2_phot)
+        energy_phot = p_phot*c
+        chi_phot = calc_chi_gamma(p0, E_f, B_f)
+        gamma_phot = np.linalg.norm(p0)/mec
+
+        print("** Case {:d} **".format(idx+1))
+        print("  initial momentum: ", p0)
+        print("  quantum parameter: {:f}".format(chi_phot))
+        print("  normalized photon energy: {:f}".format(gamma_phot))
+        print("  timestep: {:f} fs".format(sim_time*1e15))
+
+        phot_data = get_spec(all_data_end, phot_name, is_photon=True)
+        ele_data = get_spec(all_data_end, ele_name, is_photon=False)
+        pos_data = get_spec(all_data_end, pos_name, is_photon=False)
+
+        p2_ele = ele_data["px"]**2 + ele_data["py"]**2 + ele_data["pz"]**2
+        p_ele = np.sqrt(p2_ele)
+        energy_ele = np.sqrt(1.0 + p2_ele/mec**2 )*mec2
+        p2_pos = pos_data["px"]**2 + pos_data["py"]**2 + pos_data["pz"]**2
+        p_pos = np.sqrt(p2_pos)
+        energy_pos = np.sqrt(1.0 + p2_pos/mec**2 )*mec2
+
+        n_lost = check_number_of_pairs(data_set_end,
+                              phot_name, ele_name, pos_name,
+                              chi_phot, gamma_phot, sim_time,
+                              initial_particle_number)
+
+        check_weights(phot_data, ele_data, pos_data)
+
+        check_momenta(phot_data, ele_data, pos_data, p0, p_ele, p_pos)
+
+        check_energy(energy_phot, energy_ele, energy_pos)
+
+        check_energy_distrib(energy_ele, energy_pos, gamma_phot, chi_phot, n_lost, NNS[idx], idx)
+
+        check_opt_depths(phot_data, ele_data, pos_data)
+
+        print("*************\n")
+
+    test_name = filename_end[:-9] # Could also be os.path.split(os.getcwd())[1]
+    checksumAPI.evaluate_checksum(test_name, filename_end)
+
+def main():
+    check()
+
+if __name__ == "__main__":
+    main()
diff --git a/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py b/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py
deleted file mode 100755
index d4b7666bd..000000000
--- a/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py
+++ /dev/null
@@ -1,50 +0,0 @@
-#! /usr/bin/env python
-
-# Copyright 2019 Luca Fedeli, Maxence Thevenet
-#
-# This file is part of WarpX.
-#
-# License: BSD-3-Clause-LBNL
-
-import yt
-import numpy as np
-import scipy.stats as st
-import sys
-sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
-import checksumAPI
-
-# This script checks if photons initialized with Breit Wheeler process enabled
-# do actually have an exponentially distributed optical depth
-
-# Tolerance
-tolerance_rel = 1e-2
-
-def check():
-    filename = sys.argv[1]
-    data_set = yt.load(filename)
-
-    all_data = data_set.all_data()
-    res_tau = all_data["photons", 'particle_optical_depth_BW']
-
-    loc, scale = st.expon.fit(res_tau)
-
-    # loc should be very close to 0, scale should be very close to 1
-
-    error_rel = np.abs(loc - 0)
-    print("error_rel for location: " + str(error_rel))
-    print("tolerance_rel: " + str(tolerance_rel))
-    assert( error_rel < tolerance_rel )
-
-    error_rel = np.abs(scale - 1)
-    print("error_rel for scale: " + str(error_rel))
-    print("tolerance_rel: " + str(tolerance_rel))
-    assert( error_rel < tolerance_rel )
-
-    test_name = filename[:-9] # Could also be os.path.split(os.getcwd())[1]
-    checksumAPI.evaluate_checksum(test_name, filename)
-
-def main():
-    check()
-
-if __name__ == "__main__":
-    main()
diff --git a/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py b/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py
deleted file mode 100755
index 876b6eac4..000000000
--- a/Examples/Modules/qed/breit_wheeler/analysis_3d_optical_depth_evolution.py
+++ /dev/null
@@ -1,145 +0,0 @@
-#! /usr/bin/env python
-
-# Copyright 2019 Luca Fedeli, Maxence Thevenet
-#
-# This file is part of WarpX.
-#
-# License: BSD-3-Clause-LBNL
-
-# -*- coding: utf-8 -*-
-
-import yt
-import numpy as np
-import sys
-import math as m
-import scipy.special as spe
-import scipy.integrate as integ
-import scipy.stats as st
-sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
-import checksumAPI
-
-# This script checks if the optical depth of photons undergoing the
-# Breit Wheeler process behaves as expected. Four populations of photons
-# are initialized with different momenta in a background EM field. The decrease
-# of the optical depth (averaged for each population) is used to estimate the
-# Breit Wheeler cross section, which is then compared with the theoretical
-# formula. Relative discrepancy should be smaller than 2%
-#
-# References:
-# 1) R. Duclous et al 2011 Plasma Phys. Control. Fusion 53 015009
-# 2) A. Gonoskov et al. 2015 Phys. Rev. E 92, 023305
-# 3) M. Lobet. PhD thesis "Effets radiatifs et d'electrodynamique
-#    quantique dans l'interaction laser-matiere ultra-relativiste"
-#    URL: https://tel.archives-ouvertes.fr/tel-01314224
-
-
-# Tolerance
-tol = 7.e-2
-
-# EM fields
-E_f = np.array([-2433321316961438, 973328526784575, 1459992790176863])
-B_f = np.array([2857142.85714286, 4285714.28571428, 8571428.57142857])
-
-# Physical constants
-electron_mass = 9.10938356e-31
-elementary_charge = 1.6021766208e-19
-speed_of_light = 299792458
-reduced_plank = 1.054571800e-34
-vacuum_permittivity =  8.854187817e-12
-fine_structure_constant =  0.0072973525664
-classical_elec_radius = (1./4./np.pi/vacuum_permittivity)*( elementary_charge**2 / (electron_mass * speed_of_light**2))
-lambda_ref = 1.0e-6
-field_reference = 2.0 * np.pi * electron_mass*speed_of_light * speed_of_light / (elementary_charge*lambda_ref)
-schwinger_field_SI = electron_mass**2 * speed_of_light**3/(reduced_plank*elementary_charge)
-schwinger_field_norm = electron_mass*speed_of_light*lambda_ref/(2.0*reduced_plank*m.pi)
-#______________
-
-# Initial parameters
-spec_names = ["p1", "p2", "p3", "p4"]
-mec = electron_mass*speed_of_light
-p_begin = {
-    "p1": np.array([2000.0,0,0])*mec,
-    "p2": np.array([0.0,5000.0,0.0])*mec,
-    "p3": np.array([0.0,0.0,10000.0])*mec,
-    "p4": np.array([57735.02691896, 57735.02691896, 57735.02691896])*mec
-}
-initial_particle_number = 16384
-#______________
-
-def calc_chi_gamma(p, E, B):
-    p = p / electron_mass / speed_of_light
-    E = E / field_reference
-    B = B * speed_of_light / field_reference
-    gamma_phot = np.linalg.norm(p)
-    c = p/gamma_phot
-    loc_field = gamma_phot * np.linalg.norm( E - np.dot(c,E)*c + np.cross(c,B))
-    return loc_field/schwinger_field_norm
-
-#Auxiliary functions
-def X(chi_phot, chi_ele):
-    if (chi_phot > chi_ele and chi_ele != 0):
-        return np.power(chi_phot/(chi_ele*(chi_phot-chi_ele)), 2./3.)
-    else:
-        return 1.0e30
-
-def T(chi_phot):
-    coeff = 1./(np.pi * np.sqrt(3.) * chi_phot * chi_phot)
-    inner = lambda x : integ.quad(lambda s: np.sqrt(s)*spe.kv(1./3., 2./3. * s**(3./2.)), x, np.inf)[0]
-    return integ.quad(lambda chi_ele:
-                      coeff*(inner(X(chi_phot, chi_ele)) -
-                      (2.0 - chi_phot*np.power(X(chi_phot, chi_ele), 3./2.))*spe.kv(2./3., 2./3. *X(chi_phot, chi_ele)**(3./2.)) )
-                      , 0, chi_phot)[0]
-#__________________
-
-# Breit-Wheeler total cross section
-def dNBW_dt(chi_phot, energy_phot):
-    energy_phot = energy_phot/electron_mass/speed_of_light/speed_of_light
-    return ((electron_mass*(speed_of_light)**2)*fine_structure_constant/reduced_plank)*(chi_phot/energy_phot)*T(chi_phot)
-#__________________
-
-def check():
-    filename_end = sys.argv[1]
-    data_set_end = yt.load(filename_end)
-
-    sim_time = data_set_end.current_time.to_value()
-    all_data_end = data_set_end.all_data()
-
-    for name in spec_names:
-        opt_depth = all_data_end[name, 'particle_optical_depth_BW']
-
-        #check that the distribution is still exponential with scale 1 and loc 0
-        opt_depth_loc, opt_depth_scale = st.expon.fit(opt_depth)
-        exp_loc = 0.0
-        exp_scale = 1.0
-        loc_discrepancy = np.abs(opt_depth_loc-exp_loc)
-        scale_discrepancy = np.abs(opt_depth_scale-exp_scale)
-        print("tolerance_rel: " + str(tol))
-        print("species " + name)
-        print("exp distrib loc tol = " + str(tol))
-        print("exp distrib loc discrepancy = " + str(loc_discrepancy))
-        assert(loc_discrepancy < tol)
-        print("exp distrib scale tol = " + str(tol))
-        print("exp distrib scale discrepancy = " + str(scale_discrepancy/exp_scale))
-        assert(scale_discrepancy/exp_scale < tol)
-        ###
-
-        #check if number of lost photons is (n0* (1 - exp(-rate*t)) )
-        dNBW_dt_theo = dNBW_dt(
-            calc_chi_gamma(p_begin[name], E_f, B_f),
-                np.linalg.norm(p_begin[name]*speed_of_light))
-        exp_lost= initial_particle_number*(1.0 - np.exp(-dNBW_dt_theo*sim_time))
-        lost =  initial_particle_number-np.size(opt_depth)
-        discrepancy_lost = np.abs(exp_lost-lost)
-        print("lost fraction tol = " + str(tol))
-        print("lost fraction discrepancy = " + str(discrepancy_lost/exp_lost))
-        assert(discrepancy_lost/exp_lost < tol)
-        ###
-
-    test_name = filename_end[:-9] # Could also be os.path.split(os.getcwd())[1]
-    checksumAPI.evaluate_checksum(test_name, filename_end)
-
-def main():
-    check()
-
-if __name__ == "__main__":
-    main()
diff --git a/Examples/Modules/qed/breit_wheeler/inputs_2d b/Examples/Modules/qed/breit_wheeler/inputs_2d
new file mode 100644
index 000000000..2fa4690e1
--- /dev/null
+++ b/Examples/Modules/qed/breit_wheeler/inputs_2d
@@ -0,0 +1,266 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  32 32
+amr.max_grid_size = 16   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1  1  1      # Is periodic?
+geometry.prob_lo     =  -0.5e-6   -0.5e-6  # physical domain
+geometry.prob_hi     =  0.5e-6   0.5e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = energy-conserving
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = p1 p2 p3 p4 ele1 ele2 ele3 ele4 pos1 pos2 pos3 pos4 dummy_phot
+particles.photon_species = p1 p2 p3 p4 dummy_phot
+#################################
+
+p1.species_type = "photon"
+p1.injection_style = "NUniformPerCell"
+p1.profile = "constant"
+p1.num_particles_per_cell_each_dim = 32 32
+p1.density = 1e19
+p1.profile = "constant"
+p1.momentum_distribution_type = "gaussian"
+p1.ux_m = 2000.0
+##########QED####################
+p1.do_qed = 1
+p1.do_qed_breit_wheeler = 1
+p1.qed_breit_wheeler_ele_product_species = ele1
+p1.qed_breit_wheeler_pos_product_species = pos1
+#################################
+
+p2.species_type = "photon"
+p2.injection_style = "NUniformPerCell"
+p2.profile = "constant"
+p2.num_particles_per_cell_each_dim = 32 32
+p2.density = 1e19
+p2.profile = "constant"
+p2.momentum_distribution_type = "gaussian"
+p2.uy_m = 5000.0
+##########QED####################
+p2.do_qed = 1
+p2.do_qed_breit_wheeler = 1
+p2.qed_breit_wheeler_ele_product_species = ele2
+p2.qed_breit_wheeler_pos_product_species = pos2
+#################################
+
+p3.species_type = "photon"
+p3.injection_style = "NUniformPerCell"
+p3.profile = "constant"
+p3.num_particles_per_cell_each_dim = 32 32
+p3.density = 1e19
+p3.profile = "constant"
+p3.momentum_distribution_type = "gaussian"
+p3.uz_m = 10000.0
+##########QED####################
+p3.do_qed = 1
+p3.do_qed_breit_wheeler = 1
+p3.qed_breit_wheeler_ele_product_species = ele3
+p3.qed_breit_wheeler_pos_product_species = pos3
+#################################
+
+p4.species_type = "photon"
+p4.injection_style = "NUniformPerCell"
+p4.profile = "constant"
+p4.num_particles_per_cell_each_dim = 32 32
+p4.density = 1e19
+p4.profile = "constant"
+p4.momentum_distribution_type = "gaussian"
+p4.ux_m = 57735.02691896
+p4.uy_m = 57735.02691896
+p4.uz_m = 57735.02691896
+##########QED####################
+p4.do_qed = 1
+p4.do_qed_breit_wheeler = 1
+p4.qed_breit_wheeler_ele_product_species = ele4
+p4.qed_breit_wheeler_pos_product_species = pos4
+#################################
+
+### PRODUCT SPECIES ###
+ele1.species_type = "electron"
+ele1.injection_style = nuniformpercell
+ele1.num_particles_per_cell_each_dim = 0 0
+ele1.profile = constant
+ele1.density = 0.0
+ele1.momentum_distribution_type = "gaussian"
+ele1.do_not_push = 1
+ele1.do_qed = 1
+ele1.do_qed_quantum_sync = 1
+ele1.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele2.species_type = "electron"
+ele2.injection_style = nuniformpercell
+ele2.num_particles_per_cell_each_dim = 1 1
+ele2.profile = constant
+ele2.density = 0.0
+ele2.momentum_distribution_type = "gaussian"
+ele2.do_not_push = 1
+ele2.do_qed = 1
+ele2.do_qed_quantum_sync = 1
+ele2.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele3.species_type = "electron"
+ele3.injection_style = nuniformpercell
+ele3.num_particles_per_cell_each_dim = 1 1
+ele3.profile = constant
+ele3.density = 0.0
+ele3.momentum_distribution_type = "gaussian"
+ele3.do_not_push = 1
+ele3.do_qed = 1
+ele3.do_qed_quantum_sync = 1
+ele3.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele4.species_type = "electron"
+ele4.injection_style = nuniformpercell
+ele4.num_particles_per_cell_each_dim = 1 1
+ele4.profile = constant
+ele4.density = 0.0
+ele4.momentum_distribution_type = "gaussian"
+ele4.do_not_push = 1
+ele4.do_qed = 1
+ele4.do_qed_quantum_sync = 1
+ele4.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos1.species_type = "positron"
+pos1.injection_style = nuniformpercell
+pos1.num_particles_per_cell_each_dim = 0 0
+pos1.profile = constant
+pos1.density = 0.0
+pos1.momentum_distribution_type = "gaussian"
+pos1.do_not_push = 1
+pos1.do_qed = 1
+pos1.do_qed_quantum_sync = 1
+pos1.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos2.species_type = "positron"
+pos2.injection_style = nuniformpercell
+pos2.num_particles_per_cell_each_dim = 0 0
+pos2.profile = constant
+pos2.density = 0.0
+pos2.momentum_distribution_type = "gaussian"
+pos2.do_not_push = 1
+pos2.do_qed = 1
+pos2.do_qed_quantum_sync = 1
+pos2.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos3.species_type = "positron"
+pos3.injection_style = nuniformpercell
+pos3.num_particles_per_cell_each_dim = 0 0
+pos3.profile = constant
+pos3.density = 0.0
+pos3.momentum_distribution_type = "gaussian"
+pos3.do_not_push = 1
+pos3.do_qed = 1
+pos3.do_qed_quantum_sync = 1
+pos3.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos4.species_type = "positron"
+pos4.injection_style = nuniformpercell
+pos4.num_particles_per_cell_each_dim = 0 0
+pos4.profile = constant
+pos4.density = 0.0
+pos4.momentum_distribution_type = "gaussian"
+pos4.do_not_push = 1
+pos4.do_qed = 1
+pos4.do_qed_quantum_sync = 1
+pos4.qed_quantum_sync_phot_product_species = dummy_phot
+
+dummy_phot.species_type = "photon"
+dummy_phot.injection_style = nuniformpercell
+dummy_phot.num_particles_per_cell_each_dim = 0 0
+dummy_phot.profile = constant
+dummy_phot.density = 0.0
+dummy_phot.momentum_distribution_type = "gaussian"
+dummy_phot.do_qed = 0
+
+#################################
+
+##########QED TABLES####################
+qed_bw.chi_min = 0.001
+
+qed_bw.lookup_table_mode = "builtin"
+
+#qed_bw.lookup_table_mode = "generate"
+#qed_bw.tab_dndt_chi_min = 0.01
+#qed_bw.tab_dndt_chi_max = 1000.0
+#qed_bw.tab_dndt_how_many = 256
+#qed_bw.tab_pair_chi_min = 0.01
+#qed_bw.tab_pair_chi_max = 1000.0
+#qed_bw.tab_pair_chi_how_many = 256
+#qed_bw.tab_pair_frac_how_many = 256
+#qed_bw.save_table_in = "bw_table"
+
+#qed_bw.lookup_table_mode = "load"
+#qed_bw.load_table_from = "bw_table"
+
+qed_qs.chi_min = 0.001
+
+qed_qs.lookup_table_mode = "builtin"
+
+qed_qs.photon_creation_energy_threshold = 2
+
+#qed_qs.lookup_table_mode = "generate"
+#qed_qs.tab_dndt_chi_min = 0.001
+#qed_qs.tab_dndt_chi_max = 1000.0
+#qed_qs.tab_dndt_how_many = 256
+#qed_qs.tab_em_chi_min = 0.001
+#qed_qs.tab_em_frac_min = 1.0e-12
+#qed_qs.tab_em_chi_max = 1000.0
+#qed_qs.tab_em_chi_how_many = 256
+#qed_qs.tab_em_frac_how_many = 256
+#qed_qs.save_table_in = "qs_table"
+
+#qed_qs.lookup_table_mode = "load"
+#qed_qs.load_table_from = "qs_table"
+#################################
+
+### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] )
+particles.E_ext_particle_init_style = "constant"
+particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863
+####
+
+### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] )
+particles.B_ext_particle_init_style = "constant"
+particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857
+####
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.period = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex
+diag1.p1.variables = ux uy uz w
+diag1.p2.variables = ux uy uz w
+diag1.p3.variables = ux uy uz w
+diag1.p4.variables = ux uy uz w
+
+diag1.ele1.variables = ux uy uz w
+diag1.ele2.variables = ux uy uz w
+diag1.ele3.variables = ux uy uz w
+diag1.ele4.variables = ux uy uz w
+
+diag1.pos1.variables = ux uy uz w
+diag1.pos2.variables = ux uy uz w
+diag1.pos3.variables = ux uy uz w
+diag1.pos4.variables = ux uy uz w
diff --git a/Examples/Modules/qed/breit_wheeler/inputs_2d_tau_init b/Examples/Modules/qed/breit_wheeler/inputs_2d_tau_init
deleted file mode 100644
index af77c97e0..000000000
--- a/Examples/Modules/qed/breit_wheeler/inputs_2d_tau_init
+++ /dev/null
@@ -1,107 +0,0 @@
-#################################
-####### GENERAL PARAMETERS ######
-#################################
-max_step = 1
-amr.n_cell =  128 128
-amr.max_grid_size = 128   # maximum size of each AMReX box, used to decompose the domain
-amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain
-geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 0     0      # Is periodic?
-geometry.prob_lo     =  -32.e-6   -32.e-6    # physical domain
-geometry.prob_hi     =  32.e-6   32.e-6
-amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
-
-#################################
-############ NUMERICS ###########
-#################################
-algo.current_deposition = esirkepov
-algo.charge_deposition = standard
-algo.field_gathering = energy-conserving
-algo.particle_pusher = boris
-interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
-interpolation.noy = 3
-interpolation.noz = 3
-warpx.verbose = 1
-warpx.do_dive_cleaning = 0
-warpx.use_filter = 1
-warpx.cfl = 1. # if 1., the time step is set to its CFL limit
-warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
-warpx.serialize_ics = 1
-
-#################################
-############ PLASMA #############
-#################################
-particles.species_names = photons ele_bw pos_bw
-particles.photon_species = photons
-#################################
-
-photons.species_type = "photon"
-photons.injection_style = "NUniformPerCell"
-photons.profile = "constant"
-photons.xmin = -30e-6
-photons.ymin = -30e-6
-photons.zmin = -30e-6
-photons.xmax = 30e-6
-photons.ymax = 30e-6
-photons.zmax = 30e-6
-photons.num_particles_per_cell_each_dim = 2 2
-photons.density = 1e19
-photons.profile = "constant"
-photons.momentum_distribution_type = "gaussian"
-photons.ux_m = 0.0
-photons.uy_m = 0.0
-photons.uz_m = 0.0
-photons.ux_th = 100.
-photons.uy_th = 100.
-photons.uz_th = 100.
-##########QED####################
-photons.do_qed = 1
-photons.do_qed_breit_wheeler = 1
-photons.qed_breit_wheeler_ele_product_species = ele_bw
-photons.qed_breit_wheeler_pos_product_species = pos_bw
-#################################
-
-### PRODUCT SPECIES ###
-ele_bw.species_type = "electron"
-ele_bw.injection_style = nuniformpercell
-ele_bw.num_particles_per_cell_each_dim = 1 1
-ele_bw.profile = constant
-ele_bw.density = 0.0
-ele_bw.momentum_distribution_type = "gaussian"
-ele_bw.xmin = 1 ## Ugly trick to avoid electrons at T=0
-ele_bw.xmax = -1 ##  Ugly trick to avoid electrons at T=0
-ele_bw.do_qed = 0
-
-pos_bw.species_type = "positron"
-pos_bw.injection_style = nuniformpercell
-pos_bw.num_particles_per_cell_each_dim = 1 1
-pos_bw.profile = constant
-pos_bw.density = 0.0
-pos_bw.momentum_distribution_type = "gaussian"
-pos_bw.xmin = 1 ## Ugly trick to avoid positrons at T=0
-pos_bw.xmax = -1 ##  Ugly trick to avoid positrons at T=0
-pos_bw.do_qed = 0
-#################################
-
-#######QED ENGINE PARAMETERS#####
-qed_bw.lookup_table_mode = "dummy_builtin"
-
-#qed_bw.lookup_table_mode = "generate"
-#qed_bw.chi_min = 0.001
-#qed_bw.tab_dndt_chi_min = 0.1
-#qed_bw.tab_dndt_chi_max = 200
-#qed_bw.tab_dndt_how_many = 64
-#qed_bw.tab_pair_chi_min = 0.01
-#qed_bw.tab_pair_chi_max = 200
-#qed_bw.tab_pair_chi_how_many = 2
-#qed_bw.tab_pair_frac_how_many = 2
-#qed_bw.save_table_in = "bw_micro_table"
-
-#qed_bw.lookup_table_mode = "load"
-#qed_bw.load_table_from = "bw_micro_table"
-#################################
-
-# Diagnostics
-diagnostics.diags_names = diag1
-diag1.period = 1
-diag1.diag_type = Full
diff --git a/Examples/Modules/qed/breit_wheeler/inputs_3d b/Examples/Modules/qed/breit_wheeler/inputs_3d
new file mode 100644
index 000000000..af479d706
--- /dev/null
+++ b/Examples/Modules/qed/breit_wheeler/inputs_3d
@@ -0,0 +1,266 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  16 16 16
+amr.max_grid_size = 16   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1  1  1      # Is periodic?
+geometry.prob_lo     =  -0.5e-6   -0.5e-6  -0.5e-6  # physical domain
+geometry.prob_hi     =  0.5e-6   0.5e-6   0.5e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = energy-conserving
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = p1 p2 p3 p4 ele1 ele2 ele3 ele4 pos1 pos2 pos3 pos4 dummy_phot
+particles.photon_species = p1 p2 p3 p4 dummy_phot
+#################################
+
+p1.species_type = "photon"
+p1.injection_style = "NUniformPerCell"
+p1.profile = "constant"
+p1.num_particles_per_cell_each_dim = 8 8 4
+p1.density = 1e19
+p1.profile = "constant"
+p1.momentum_distribution_type = "gaussian"
+p1.ux_m = 2000.0
+##########QED####################
+p1.do_qed = 1
+p1.do_qed_breit_wheeler = 1
+p1.qed_breit_wheeler_ele_product_species = ele1
+p1.qed_breit_wheeler_pos_product_species = pos1
+#################################
+
+p2.species_type = "photon"
+p2.injection_style = "NUniformPerCell"
+p2.profile = "constant"
+p2.num_particles_per_cell_each_dim = 8 8 4
+p2.density = 1e19
+p2.profile = "constant"
+p2.momentum_distribution_type = "gaussian"
+p2.uy_m = 5000.0
+##########QED####################
+p2.do_qed = 1
+p2.do_qed_breit_wheeler = 1
+p2.qed_breit_wheeler_ele_product_species = ele2
+p2.qed_breit_wheeler_pos_product_species = pos2
+#################################
+
+p3.species_type = "photon"
+p3.injection_style = "NUniformPerCell"
+p3.profile = "constant"
+p3.num_particles_per_cell_each_dim = 8 8 4
+p3.density = 1e19
+p3.profile = "constant"
+p3.momentum_distribution_type = "gaussian"
+p3.uz_m = 10000.0
+##########QED####################
+p3.do_qed = 1
+p3.do_qed_breit_wheeler = 1
+p3.qed_breit_wheeler_ele_product_species = ele3
+p3.qed_breit_wheeler_pos_product_species = pos3
+#################################
+
+p4.species_type = "photon"
+p4.injection_style = "NUniformPerCell"
+p4.profile = "constant"
+p4.num_particles_per_cell_each_dim = 8 8 4
+p4.density = 1e19
+p4.profile = "constant"
+p4.momentum_distribution_type = "gaussian"
+p4.ux_m = 57735.02691896
+p4.uy_m = 57735.02691896
+p4.uz_m = 57735.02691896
+##########QED####################
+p4.do_qed = 1
+p4.do_qed_breit_wheeler = 1
+p4.qed_breit_wheeler_ele_product_species = ele4
+p4.qed_breit_wheeler_pos_product_species = pos4
+#################################
+
+### PRODUCT SPECIES ###
+ele1.species_type = "electron"
+ele1.injection_style = nuniformpercell
+ele1.num_particles_per_cell_each_dim = 0 0
+ele1.profile = constant
+ele1.density = 0.0
+ele1.momentum_distribution_type = "gaussian"
+ele1.do_not_push = 1
+ele1.do_qed = 1
+ele1.do_qed_quantum_sync = 1
+ele1.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele2.species_type = "electron"
+ele2.injection_style = nuniformpercell
+ele2.num_particles_per_cell_each_dim = 0 0
+ele2.profile = constant
+ele2.density = 0.0
+ele2.momentum_distribution_type = "gaussian"
+ele2.do_not_push = 1
+ele2.do_qed = 1
+ele2.do_qed_quantum_sync = 1
+ele2.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele3.species_type = "electron"
+ele3.injection_style = nuniformpercell
+ele3.num_particles_per_cell_each_dim = 0 0
+ele3.profile = constant
+ele3.density = 0.0
+ele3.momentum_distribution_type = "gaussian"
+ele3.do_not_push = 1
+ele3.do_qed = 1
+ele3.do_qed_quantum_sync = 1
+ele3.qed_quantum_sync_phot_product_species = dummy_phot
+
+ele4.species_type = "electron"
+ele4.injection_style = nuniformpercell
+ele4.num_particles_per_cell_each_dim = 0 0
+ele4.profile = constant
+ele4.density = 0.0
+ele4.momentum_distribution_type = "gaussian"
+ele4.do_not_push = 1
+ele4.do_qed = 1
+ele4.do_qed_quantum_sync = 1
+ele4.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos1.species_type = "positron"
+pos1.injection_style = nuniformpercell
+pos1.num_particles_per_cell_each_dim = 0 0
+pos1.profile = constant
+pos1.density = 0.0
+pos1.momentum_distribution_type = "gaussian"
+pos1.do_not_push = 1
+pos1.do_qed = 1
+pos1.do_qed_quantum_sync = 1
+pos1.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos2.species_type = "positron"
+pos2.injection_style = nuniformpercell
+pos2.num_particles_per_cell_each_dim = 0 0
+pos2.profile = constant
+pos2.density = 0.0
+pos2.momentum_distribution_type = "gaussian"
+pos2.do_not_push = 1
+pos2.do_qed = 1
+pos2.do_qed_quantum_sync = 1
+pos2.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos3.species_type = "positron"
+pos3.injection_style = nuniformpercell
+pos3.num_particles_per_cell_each_dim = 0 0
+pos3.profile = constant
+pos3.density = 0.0
+pos3.momentum_distribution_type = "gaussian"
+pos3.do_not_push = 1
+pos3.do_qed = 1
+pos3.do_qed_quantum_sync = 1
+pos3.qed_quantum_sync_phot_product_species = dummy_phot
+
+pos4.species_type = "positron"
+pos4.injection_style = nuniformpercell
+pos4.num_particles_per_cell_each_dim = 0 0
+pos4.profile = constant
+pos4.density = 0.0
+pos4.momentum_distribution_type = "gaussian"
+pos4.do_not_push = 1
+pos4.do_qed = 1
+pos4.do_qed_quantum_sync = 1
+pos4.qed_quantum_sync_phot_product_species = dummy_phot
+
+dummy_phot.species_type = "photon"
+dummy_phot.injection_style = nuniformpercell
+dummy_phot.num_particles_per_cell_each_dim = 0 0
+dummy_phot.profile = constant
+dummy_phot.density = 0.0
+dummy_phot.momentum_distribution_type = "gaussian"
+dummy_phot.do_qed = 0
+
+#################################
+
+##########QED TABLES####################
+qed_bw.chi_min = 0.001
+
+qed_bw.lookup_table_mode = "builtin"
+
+#qed_bw.lookup_table_mode = "generate"
+#qed_bw.tab_dndt_chi_min = 0.01
+#qed_bw.tab_dndt_chi_max = 1000.0
+#qed_bw.tab_dndt_how_many = 256
+#qed_bw.tab_pair_chi_min = 0.01
+#qed_bw.tab_pair_chi_max = 1000.0
+#qed_bw.tab_pair_chi_how_many = 256
+#qed_bw.tab_pair_frac_how_many = 256
+#qed_bw.save_table_in = "bw_table"
+
+#qed_bw.lookup_table_mode = "load"
+#qed_bw.load_table_from = "bw_table"
+
+qed_qs.chi_min = 0.001
+
+qed_qs.lookup_table_mode = "builtin"
+
+qed_qs.photon_creation_energy_threshold = 2
+
+#qed_qs.lookup_table_mode = "generate"
+#qed_qs.tab_dndt_chi_min = 0.001
+#qed_qs.tab_dndt_chi_max = 1000.0
+#qed_qs.tab_dndt_how_many = 256
+#qed_qs.tab_em_chi_min = 0.001
+#qed_qs.tab_em_frac_min = 1.0e-12
+#qed_qs.tab_em_chi_max = 1000.0
+#qed_qs.tab_em_chi_how_many = 256
+#qed_qs.tab_em_frac_how_many = 256
+#qed_qs.save_table_in = "qs_table"
+
+#qed_qs.lookup_table_mode = "load"
+#qed_qs.load_table_from = "qs_table"
+#################################
+
+### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] )
+particles.E_ext_particle_init_style = "constant"
+particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863
+####
+
+### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] )
+particles.B_ext_particle_init_style = "constant"
+particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857
+####
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.period = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex
+diag1.p1.variables = ux uy uz w
+diag1.p2.variables = ux uy uz w
+diag1.p3.variables = ux uy uz w
+diag1.p4.variables = ux uy uz w
+
+diag1.ele1.variables = ux uy uz w
+diag1.ele2.variables = ux uy uz w
+diag1.ele3.variables = ux uy uz w
+diag1.ele4.variables = ux uy uz w
+
+diag1.pos1.variables = ux uy uz w
+diag1.pos2.variables = ux uy uz w
+diag1.pos3.variables = ux uy uz w
+diag1.pos4.variables = ux uy uz w
diff --git a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution b/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution
deleted file mode 100644
index db4ec895e..000000000
--- a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution
+++ /dev/null
@@ -1,200 +0,0 @@
-#################################
-####### GENERAL PARAMETERS ######
-#################################
-max_step = 5
-amr.n_cell =  32 32 16
-amr.max_grid_size = 32   # maximum size of each AMReX box, used to decompose the domain
-amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
-geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 1  1  1      # Is periodic?
-geometry.prob_lo     =  -0.5e-6   -0.5e-6  -0.25e-6  # physical domain
-geometry.prob_hi     =  0.5e-6   0.5e-6   0.25e-6
-amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
-
-#################################
-############ NUMERICS ###########
-#################################
-algo.current_deposition = esirkepov
-algo.charge_deposition = standard
-algo.field_gathering = energy-conserving
-algo.particle_pusher = boris
-interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
-interpolation.noy = 3
-interpolation.noz = 3
-warpx.verbose = 1
-warpx.do_dive_cleaning = 0
-warpx.use_filter = 1
-warpx.cfl = 1. # if 1., the time step is set to its CFL limit
-warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
-warpx.serialize_ics = 1
-
-#################################
-############ PLASMA #############
-#################################
-particles.species_names = p1 p2 p3 p4 ele_bw pos_bw
-particles.photon_species = p1 p2 p3 p4
-#################################
-
-p1.species_type = "photon"
-p1.injection_style = "NUniformPerCell"
-p1.profile = "constant"
-p1.xmin = -0.6e-6
-p1.ymin = -0.6e-6
-p1.zmin = -0.6e-6
-p1.xmax = 0.6e6
-p1.ymax = 0.6e6
-p1.zmax = 0.6e6
-p1.num_particles_per_cell_each_dim = 1 1 1
-p1.density = 1e19
-p1.profile = "constant"
-p1.momentum_distribution_type = "gaussian"
-p1.ux_m = 2000.0
-p1.uy_m = 0.0
-p1.uz_m = 0.0
-p1.ux_th = 0.
-p1.uy_th = 0.
-p1.uz_th = 0.
-##########QED####################
-p1.do_qed = 1
-p1.do_qed_breit_wheeler = 1
-p1.qed_breit_wheeler_ele_product_species = ele_bw
-p1.qed_breit_wheeler_pos_product_species = pos_bw
-#################################
-
-p2.species_type = "photon"
-p2.injection_style = "NUniformPerCell"
-p2.profile = "constant"
-p2.xmin = -0.6e-6
-p2.ymin = -0.6e-6
-p2.zmin = -0.6e-6
-p2.xmax = 0.6e6
-p2.ymax = 0.6e6
-p2.zmax = 0.6e6
-p2.num_particles_per_cell_each_dim = 1 1 1
-p2.density = 1e19
-p2.profile = "constant"
-p2.momentum_distribution_type = "gaussian"
-p2.ux_m = 0.0
-p2.uy_m = 5000.0
-p2.uz_m = 0.0
-p2.ux_th = 0.
-p2.uy_th = 0.
-p2.uz_th = 0.
-##########QED####################
-p2.do_qed = 1
-p2.do_qed_breit_wheeler = 1
-p2.qed_breit_wheeler_ele_product_species = ele_bw
-p2.qed_breit_wheeler_pos_product_species = pos_bw
-#################################
-
-p3.species_type = "photon"
-p3.injection_style = "NUniformPerCell"
-p3.profile = "constant"
-p3.xmin = -0.6e-6
-p3.ymin = -0.6e-6
-p3.zmin = -0.6e-6
-p3.xmax = 0.6e6
-p3.ymax = 0.6e6
-p3.zmax = 0.6e6
-p3.num_particles_per_cell_each_dim = 1 1 1
-p3.density = 1e19
-p3.profile = "constant"
-p3.momentum_distribution_type = "gaussian"
-p3.ux_m = 0.0
-p3.uy_m = 0.0
-p3.uz_m = 10000.0
-p3.ux_th = 0.
-p3.uy_th = 0.
-p3.uz_th = 0.
-##########QED####################
-p3.do_qed = 1
-p3.do_qed_breit_wheeler = 1
-p3.qed_breit_wheeler_ele_product_species = ele_bw
-p3.qed_breit_wheeler_pos_product_species = pos_bw
-#################################
-
-p4.species_type = "photon"
-p4.injection_style = "NUniformPerCell"
-p4.profile = "constant"
-p4.xmin = -0.6e-6
-p4.ymin = -0.6e-6
-p4.zmin = -0.6e-6
-p4.xmax = 0.6e6
-p4.ymax = 0.6e6
-p4.zmax = 0.6e6
-p4.num_particles_per_cell_each_dim = 1 1 1
-p4.density = 1e19
-p4.profile = "constant"
-p4.momentum_distribution_type = "gaussian"
-p4.ux_m = 57735.02691896
-p4.uy_m = 57735.02691896
-p4.uz_m = 57735.02691896
-p4.ux_th = 0.
-p4.uy_th = 0.
-p4.uz_th = 0.
-##########QED####################
-p4.do_qed = 1
-p4.do_qed_breit_wheeler = 1
-p4.qed_breit_wheeler_ele_product_species = ele_bw
-p4.qed_breit_wheeler_pos_product_species = pos_bw
-#################################
-
-### PRODUCT SPECIES ###
-ele_bw.species_type = "electron"
-ele_bw.injection_style = nuniformpercell
-ele_bw.num_particles_per_cell_each_dim = 1 1
-ele_bw.profile = constant
-ele_bw.density = 0.0
-ele_bw.momentum_distribution_type = "gaussian"
-ele_bw.xmin = 1 ## Ugly trick to avoid electrons at T=0
-ele_bw.xmax = -1 ##  Ugly trick to avoid electrons at T=0
-ele_bw.do_qed = 0
-
-pos_bw.species_type = "positron"
-pos_bw.injection_style = nuniformpercell
-pos_bw.num_particles_per_cell_each_dim = 1 1
-pos_bw.profile = constant
-pos_bw.density = 0.0
-pos_bw.momentum_distribution_type = "gaussian"
-pos_bw.xmin = 1 ## Ugly trick to avoid positrons at T=0
-pos_bw.xmax = -1 ##  Ugly trick to avoid positrons at T=0
-pos_bw.do_qed = 0
-#################################
-
-##########QED TABLES####################
-qed_bw.lookup_table_mode = "dummy_builtin"
-
-#qed_bw.lookup_table_mode = "generate"
-#qed_bw.chi_min = 0.001
-#qed_bw.tab_dndt_chi_min = 0.1
-#qed_bw.tab_dndt_chi_max = 200
-#qed_bw.tab_dndt_how_many = 64
-#qed_bw.tab_pair_chi_min = 0.01
-#qed_bw.tab_pair_chi_max = 200
-#qed_bw.tab_pair_chi_how_many = 2
-#qed_bw.tab_pair_frac_how_many = 2
-#qed_bw.save_table_in = "bw_micro_table"
-
-#qed_bw.lookup_table_mode = "load"
-#qed_bw.load_table_from = "bw_micro_table"
-#################################
-
-### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] )
-particles.E_ext_particle_init_style = "constant"
-particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863
-####
-
-### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] )
-particles.B_ext_particle_init_style = "constant"
-particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857
-####
-
-# Diagnostics
-diagnostics.diags_names = diag1
-diag1.period = 1
-diag1.diag_type = Full
-diag1.fields_to_plot = Ex
-diag1.p1.variables = ux uy uz
-diag1.p2.variables = ux uy uz
-diag1.p3.variables = ux uy uz
-diag1.p4.variables = ux uy uz
diff --git a/Examples/Modules/qed/quantum_synchrotron/analysis.py b/Examples/Modules/qed/quantum_synchrotron/analysis.py
new file mode 100755
index 000000000..5c2f778a1
--- /dev/null
+++ b/Examples/Modules/qed/quantum_synchrotron/analysis.py
@@ -0,0 +1,298 @@
+#! /usr/bin/env python
+
+# Copyright 2019 Luca Fedeli, Maxence Thevenet
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# -*- coding: utf-8 -*-
+
+import yt
+import numpy as np
+import sys
+import scipy.special as spe
+import scipy.integrate as integ
+import scipy.stats as st
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+import matplotlib.pyplot as plt
+
+# This script performs detailed checks of the Quantum Synchrotron photon emission process.
+# Two electron populations and two positron populations are initialized with different momenta in different
+# directions in a background EM field (with non-zero components along each direction).
+# Specifically the script checks that:
+#
+# - The expected number of generated photons n_phot is in agreement with theory.
+#   The maximum tolerated error is 5*sqrt(n_phot), except for the last 8 points,
+#   for which a higher tolerance is used (this is due to the fact that the resolution
+#   of the builtin table is quite limited).
+# - The weight of the generated particles is equal to the weight of the photon
+# - The generated particles are emitted in the right direction
+# - The energy distribution of the generated particles is in agreement with theory
+# - The optical depths of the product species are correctly initialized (QED effects are
+#   enabled for product species too).
+#
+# More details on the theoretical formulas used in this script can be found in
+# the Jupyter notebook picsar/src/multi_physics/QED_tests/validation/validation.ipynb
+#
+# References:
+# 1) R. Duclous et al 2011 Plasma Phys. Control. Fusion 53 015009
+# 2) A. Gonoskov et al. 2015 Phys. Rev. E 92, 023305
+# 3) M. Lobet. PhD thesis "Effets radiatifs et d'electrodynamique
+#    quantique dans l'interaction laser-matiere ultra-relativiste"
+#    URL: https://tel.archives-ouvertes.fr/tel-01314224
+
+
+# Tolerances
+tol = 1.e-8
+tol_red = 1.e-2
+
+# Physical constants (from CODATA 2018, see: https://physics.nist.gov/cuu/Constants/index.html )
+me = 9.1093837015e-31 #electron mass
+c = 299792458 #speed of light
+hbar = 6.62607015e-34/(2*np.pi) #reduced Plank constant
+fine_structure = 7.2973525693e-3 #fine structure constant
+qe = 1.602176634e-19#elementary charge
+E_s = (me**2 * c**3)/(qe * hbar) #Schwinger E field
+B_s = E_s/c #Schwinger B field
+
+mec = me*c
+mec2 = mec*c
+#______________
+
+# Initial parameters
+spec_names = ["p1", "p2", "p3", "p4"]
+spec_names_phot = ["qsp_1", "qsp_2", "qsp_3", "qsp_4"]
+initial_momenta = [
+    np.array([10.0,0,0])*mec,
+    np.array([0.0,100.0,0.0])*mec,
+    np.array([0.0,0.0,1000.0])*mec,
+    np.array([5773.502691896, 5773.502691896, 5773.502691896])*mec
+]
+csign = [-1,-1,1,1]
+initial_particle_number = 1048576
+
+E_f = np.array([-2433321316961438., 973328526784575., 1459992790176863.])
+B_f = np.array([2857142.85714286, 4285714.28571428, 8571428.57142857])
+
+NNS = [64,64,64,64] #bins for energy distribution comparison.
+#______________
+
+def calc_chi_part(p, E, B):
+    gamma_part = np.sqrt(1.0 + np.dot(p,p)/mec**2)
+    v = p/(gamma_part*me)
+    EpvvecB = E + np.cross(v,B)
+    vdotEoverc = np.dot(v,E)/c
+    ff = np.sqrt(np.dot(EpvvecB,EpvvecB) - np.dot(vdotEoverc,vdotEoverc))
+    return gamma_part*ff/E_s
+
+#Auxiliary functions
+@np.vectorize
+def IC_inner_alternative(y):
+    ff = lambda x : np.exp(-y*(1+(4*x**2)/3)*np.sqrt(1+x*x/3))*(9+36*x**2 + 16*x**4)/(3 + 4*x**2)/np.sqrt(1+(x**2)/3)
+    return integ.quad(ff, 0, np.inf)[0]/np.sqrt(3)
+
+def IC_Y(chi_ele, xi):
+    div = (chi_ele*(1-xi))
+    div = np.where(np.logical_and(xi < 1, chi_ele != 0), div, 1.0);
+    res = (2/3)*np.where(np.logical_and(xi < 1, chi_ele != 0), xi/div, np.inf)
+    return res
+
+def IC_S(chi_ele, xi):
+    Y = IC_Y(chi_ele, xi)
+    coeff = np.sqrt(3)/2.0/np.pi
+    first = IC_inner_alternative(Y)
+    div = np.where(xi == 1, 1.0, 1.0/(1-xi)  )
+    res = np.where(np.logical_or(xi == 1, xi == 0), 0.0,
+        coeff*xi*( first  + (xi**2 * spe.kv(2./3.,Y)*div )  ) )
+    return res
+
+def IC_SXI(chi_ele, xi):
+    div = np.where(xi != 0, xi, 1.0)
+    return np.where(xi != 0, IC_S(chi_ele, xi)/div, np.inf)
+
+@np.vectorize
+def IC_G(chi_ele):
+    return integ.quad(lambda xi: IC_SXI(chi_ele, xi), 0, 1)[0]
+
+def small_diff(vv, val):
+    if(val != 0.0):
+        return np.max(np.abs((vv - val)/val)) < tol
+    else:
+        return np.max(np.abs(vv)) < tol
+
+def boris(pp, dt, charge_sign):
+    econst = 0.5*qe*dt*charge_sign/me
+    u = pp/(me)
+    u += econst*E_f
+    inv_gamma = 1/np.sqrt(1 + np.dot(u,u)/c**2)
+    t = econst*B_f*inv_gamma
+    s = 2*t/(1 + np.dot(t,t))
+    u_p = u + np.cross(u,t)
+    u += np.cross(u_p, s)
+    u += econst*E_f
+    return u *me
+#__________________
+
+# Quantum Synchrotron total and differential cross sections
+def QS_dN_dt(chi_ele, gamma_ele):
+    coeff_IC = (2./3.) * fine_structure * me*c**2/hbar
+    return coeff_IC*IC_G(chi_ele)/gamma_ele
+
+def QS_d2N_dt_dchi(chi, gamma_ele, chi_phot):
+    coeff_IC = (2./3.) * fine_structure * me*c**2/hbar
+    return coeff_IC*IC_S(chi, chi_phot/chi)/chi_phot/gamma_ele
+#__________________
+
+# Get data for a species
+def get_spec(ytdata, specname, is_photon):
+    px = ytdata[specname,"particle_momentum_x"].v
+    pz = ytdata[specname,"particle_momentum_z"].v
+    py = ytdata[specname,"particle_momentum_y"].v
+
+    w = ytdata[specname,"particle_weighting"].v
+
+    if (is_photon):
+        opt = ytdata[specname,"particle_optical_depth_BW"].v
+    else:
+        opt = ytdata[specname,"particle_optical_depth_QSR"].v
+
+    return {"px" : px, "py" : py, "pz" : pz, "w" : w, "opt" : opt}
+
+# Individual tests
+def check_number_of_photons(ytdataset, part_name, phot_name, chi_part, gamma_part, dt, particle_number):
+    dNQS_dt_theo = QS_dN_dt(chi_part, gamma_part)
+    expected_photons = (1.-np.exp(-dNQS_dt_theo*dt))*particle_number
+    expected_photons_tolerance = 5.0*np.sqrt(expected_photons)
+    n_phot = ytdataset.particle_type_counts[phot_name]
+    assert( np.abs(n_phot-expected_photons) < expected_photons_tolerance)
+    print("  [OK] generated photons number is within expectations")
+    return n_phot
+
+def check_weights(part_data, phot_data):
+    assert(np.all(part_data["w"] == part_data["w"][0]))
+    assert(np.all(phot_data["w"]  == part_data["w"][0]))
+    print("  [OK] particles weights are the expected ones")
+
+def check_momenta(phot_data, p_phot, p0):
+    pdir = p0/np.linalg.norm(p0)
+    assert(small_diff(phot_data["px"]/p_phot, pdir[0]))
+    assert(small_diff(phot_data["py"]/p_phot, pdir[1]))
+    assert(small_diff(phot_data["pz"]/p_phot, pdir[2]))
+    print("  [OK] photons move along the initial particle direction")
+
+def check_opt_depths(part_data, phot_data):
+    data = (part_data, phot_data)
+    for dd in data:
+        loc, scale = st.expon.fit(dd["opt"])
+        assert( np.abs(loc - 0) < tol_red )
+        assert( np.abs(scale - 1) < tol_red )
+    print("  [OK] optical depth distributions are still exponential")
+
+def check_energy_distrib(gamma_phot, chi_part,
+        gamma_part, n_phot, NN, idx, do_plot=False):
+    gamma_phot_min = 1e-12*gamma_part
+    gamma_phot_max = gamma_part
+    h_log_gamma_phot, c_gamma_phot = np.histogram(gamma_phot, bins=np.logspace(np.log10(gamma_phot_min),np.log10(gamma_phot_max),NN+1))
+
+    cchi_phot_min = chi_part*(gamma_phot_min)/(gamma_part-1)
+    cchi_phot_max = chi_part*(gamma_phot_max)/(gamma_part-1)
+
+    #Rudimentary integration over npoints for each bin
+    npoints= 20
+    aux_chi = np.logspace(np.log10(cchi_phot_min),np.log10(cchi_phot_max), NN*npoints)
+    distrib = QS_d2N_dt_dchi(chi_part, gamma_part, aux_chi)*aux_chi
+    distrib = np.sum(distrib.reshape(-1, npoints),1)
+    distrib = n_phot*distrib/np.sum(distrib)
+
+    if do_plot :
+        # Visual comparison of distributions
+        c_gamma_phot = np.exp(0.5*(np.log(c_gamma_phot[1:])+np.log(c_gamma_phot[:-1])))
+        plt.clf()
+
+        fig, (ax1, ax2) = plt.subplots(1, 2)
+        fig.suptitle("χ_particle = {:f}".format(chi_part))
+        ax1.plot(c_gamma_phot, distrib,label="theory")
+        ax1.loglog(c_gamma_phot, h_log_gamma_phot,label="QSR photons")
+        ax1.set_xlim(1e-12*(gamma_part-1),gamma_part-1)
+        ax1.set_ylim(1,1e5)
+        ax2.plot(c_gamma_phot, distrib,label="theory")
+        ax2.semilogy(c_gamma_phot, h_log_gamma_phot,label="QSR photons")
+        ax2.set_ylim(1,1e5)
+        ax2.set_xlim(1e-12*(gamma_part-1),gamma_part-1)
+        ax1.set_xlabel("γ_photon")
+        ax1.set_xlabel("N")
+        ax2.set_xlabel("γ_photon")
+        ax2.set_xlabel("N")
+        plt.legend()
+        plt.savefig("case_{:d}".format(idx+1))
+
+    discr = np.abs(h_log_gamma_phot-distrib)
+
+    max_discr = np.sqrt(distrib)*5.0
+    # Use a higer tolerance for the last 8 points (this is due to limitations
+    # of the builtin table)
+    max_discr[-8:] *= 2.0
+    assert(np.all( discr < max_discr ))
+
+    print("  [OK] energy distribution is within expectations")
+
+#__________________
+
+def check():
+    filename_end = sys.argv[1]
+    data_set_end = yt.load(filename_end)
+
+    sim_time = data_set_end.current_time.to_value()
+    all_data_end = data_set_end.all_data()
+
+    for idx in range(4):
+        part_name = spec_names[idx]
+        phot_name  = spec_names_phot[idx]
+        t_pi        = initial_momenta[idx]
+        pm = boris(t_pi,-sim_time*0.5,csign[idx])
+        p0 = boris(pm,sim_time*1.0,csign[idx])
+
+        p2_part = p0[0]**2 + p0[1]**2 + p0[2]**2
+        energy_part = np.sqrt(mec2**2 + p2_part*c**2)
+        gamma_part = energy_part/mec2
+        chi_part = calc_chi_part(p0, E_f, B_f)
+
+        print("** Case {:d} **".format(idx+1))
+        print("  initial momentum: ", t_pi)
+        print("  quantum parameter: {:f}".format(chi_part))
+        print("  normalized particle energy: {:f}".format(gamma_part))
+        print("  timestep: {:f} fs".format(sim_time*1e15))
+
+        part_data_final = get_spec(all_data_end, part_name, is_photon=False)
+        phot_data = get_spec(all_data_end, phot_name, is_photon=True)
+
+        p_phot = np.sqrt(phot_data["px"]**2 + phot_data["py"]**2 + phot_data["pz"]**2)
+        energy_phot = p_phot*c
+        gamma_phot = energy_phot/mec2
+
+        n_phot = check_number_of_photons(data_set_end,
+                              part_name, phot_name,
+                              chi_part, gamma_part, sim_time,
+                              initial_particle_number)
+
+        check_weights(part_data_final, phot_data)
+
+        check_momenta(phot_data, p_phot, p0)
+
+        check_energy_distrib(gamma_phot, chi_part, gamma_part, n_phot, NNS[idx], idx)
+
+        check_opt_depths(part_data_final, phot_data)
+
+        print("*************\n")
+
+    test_name = filename_end[:-9] # Could also be os.path.split(os.getcwd())[1]
+    checksumAPI.evaluate_checksum(test_name, filename_end)
+
+def main():
+    check()
+
+if __name__ == "__main__":
+    main()
diff --git a/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py b/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py
deleted file mode 100755
index b21441fb7..000000000
--- a/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py
+++ /dev/null
@@ -1,60 +0,0 @@
-#! /usr/bin/env python
-
-# Copyright 2019 Luca Fedeli, Maxence Thevenet
-#
-# This file is part of WarpX.
-#
-# License: BSD-3-Clause-LBNL
-
-import yt
-import numpy as np
-import scipy.stats as st
-import sys
-sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
-import checksumAPI
-
-# This script checks if electrons and positrons initialized with
-# Quantum Synchrotron process enabled
-# do actually have an exponentially distributed optical depth
-
-# Tolerance
-tolerance_rel = 1e-2
-print("tolerance_rel: " + str(tolerance_rel))
-
-def check():
-    filename = sys.argv[1]
-    data_set = yt.load(filename)
-
-    all_data = data_set.all_data()
-    res_ele_tau = all_data["electrons", 'particle_optical_depth_QSR']
-    res_pos_tau = all_data["positrons", 'particle_optical_depth_QSR']
-
-    loc_ele, scale_ele = st.expon.fit(res_ele_tau)
-    loc_pos, scale_pos = st.expon.fit(res_pos_tau)
-
-    # loc should be very close to 0, scale should be very close to 1
-    error_rel = np.abs(loc_ele - 0)
-    print("error_rel loc_ele: " + str(error_rel))
-    assert( error_rel < tolerance_rel )
-
-    error_rel = np.abs(loc_pos - 0)
-    print("error_rel loc_pos: " + str(error_rel))
-    assert( error_rel < tolerance_rel )
-
-    error_rel = np.abs(scale_ele - 1)
-    print("error_rel scale_ele: " + str(error_rel))
-    assert( error_rel < tolerance_rel )
-
-    error_rel = np.abs(scale_pos - 1)
-    print("error_rel scale_pos: " + str(error_rel))
-    assert( error_rel < tolerance_rel )
-
-    test_name = filename[:-9] # Could also be os.path.split(os.getcwd())[1]
-    checksumAPI.evaluate_checksum(test_name, filename)
-
-def main():
-    check()
-
-if __name__ == "__main__":
-    main()
-
diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs_2d b/Examples/Modules/qed/quantum_synchrotron/inputs_2d
new file mode 100644
index 000000000..c57e04793
--- /dev/null
+++ b/Examples/Modules/qed/quantum_synchrotron/inputs_2d
@@ -0,0 +1,224 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  32 32
+amr.max_grid_size = 16   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1  1  1      # Is periodic?
+geometry.prob_lo     =  -0.25e-6   -0.25e-6  # physical domain
+geometry.prob_hi     =  0.25e-6   0.25e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = energy-conserving
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = p1 p2 p3 p4 qsp_1 qsp_2 qsp_3 qsp_4 dummy_ele dummy_pos
+particles.photon_species = qsp_1 qsp_2 qsp_3 qsp_4
+#################################
+
+p1.species_type = "electron"
+p1.injection_style = "NUniformPerCell"
+p1.profile = "constant"
+p1.num_particles_per_cell_each_dim = 32 32
+p1.density = 1e1
+p1.profile = "constant"
+p1.momentum_distribution_type = "gaussian"
+p1.ux_m = 10.0
+##########QED####################
+p1.do_qed = 1
+p1.do_qed_quantum_sync = 1
+p1.qed_quantum_sync_phot_product_species = qsp_1
+#################################
+
+p2.species_type = "electron"
+p2.injection_style = "NUniformPerCell"
+p2.profile = "constant"
+p2.num_particles_per_cell_each_dim = 32 32
+p2.density = 1e1
+p2.profile = "constant"
+p2.momentum_distribution_type = "gaussian"
+p2.uy_m = 100.0
+##########QED####################
+p2.do_qed = 1
+p2.do_qed_quantum_sync = 1
+p2.qed_quantum_sync_phot_product_species = qsp_2
+#################################
+
+p3.species_type = "positron"
+p3.injection_style = "NUniformPerCell"
+p3.profile = "constant"
+p3.num_particles_per_cell_each_dim = 32 32
+p3.density = 1e1
+p3.profile = "constant"
+p3.momentum_distribution_type = "gaussian"
+p3.uz_m = 1000.0
+##########QED####################
+p3.do_qed = 1
+p3.do_qed_quantum_sync = 1
+p3.qed_quantum_sync_phot_product_species = qsp_3
+#################################
+
+p4.species_type = "positron"
+p4.injection_style = "NUniformPerCell"
+p4.profile = "constant"
+p4.num_particles_per_cell_each_dim = 32 32
+p4.density = 1e1
+p4.profile = "constant"
+p4.momentum_distribution_type = "gaussian"
+p4.ux_m = 5773.502691896
+p4.uy_m = 5773.502691896
+p4.uz_m = 5773.502691896
+##########QED####################
+p4.do_qed = 1
+p4.do_qed_quantum_sync = 1
+p4.qed_quantum_sync_phot_product_species = qsp_4
+#################################
+
+### PRODUCT SPECIES ###
+qsp_1.species_type = "photon"
+qsp_1.injection_style = nuniformpercell
+qsp_1.num_particles_per_cell_each_dim = 0 0
+qsp_1.profile = constant
+qsp_1.density = 0.0
+qsp_1.momentum_distribution_type = "gaussian"
+qsp_1.do_qed = 1
+qsp_1.do_qed_breit_wheeler = 1
+qsp_1.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_1.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_2.species_type = "photon"
+qsp_2.injection_style = nuniformpercell
+qsp_2.num_particles_per_cell_each_dim = 0 0
+qsp_2.profile = constant
+qsp_2.density = 0.0
+qsp_2.momentum_distribution_type = "gaussian"
+qsp_2.do_qed = 1
+qsp_2.do_qed_breit_wheeler = 1
+qsp_2.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_2.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_3.species_type = "photon"
+qsp_3.injection_style = nuniformpercell
+qsp_3.num_particles_per_cell_each_dim = 0 0
+qsp_3.profile = constant
+qsp_3.density = 0.0
+qsp_3.momentum_distribution_type = "gaussian"
+qsp_3.do_qed = 1
+qsp_3.do_qed_breit_wheeler = 1
+qsp_3.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_3.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_4.species_type = "photon"
+qsp_4.injection_style = nuniformpercell
+qsp_4.num_particles_per_cell_each_dim = 0 0
+qsp_4.profile = constant
+qsp_4.density = 0.0
+qsp_4.momentum_distribution_type = "gaussian"
+qsp_4.do_qed = 1
+qsp_4.do_qed_breit_wheeler = 1
+qsp_4.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_4.qed_breit_wheeler_pos_product_species = dummy_pos
+
+#################################
+
+dummy_ele.species_type = "electron"
+dummy_ele.injection_style = nuniformpercell
+dummy_ele.num_particles_per_cell_each_dim = 0 0
+dummy_ele.profile = constant
+dummy_ele.density = 0.0
+dummy_ele.momentum_distribution_type = "gaussian"
+dummy_ele.do_qed = 0
+
+dummy_pos.species_type = "positron"
+dummy_pos.injection_style = nuniformpercell
+dummy_pos.num_particles_per_cell_each_dim = 0 0
+dummy_pos.profile = constant
+dummy_pos.density = 0.0
+dummy_pos.momentum_distribution_type = "gaussian"
+dummy_pos.do_qed = 0
+
+
+#################################
+
+##########QED TABLES####################
+qed_bw.chi_min = 0.001
+
+qed_bw.lookup_table_mode = "builtin"
+
+#qed_bw.lookup_table_mode = "generate"
+#qed_bw.tab_dndt_chi_min = 0.01
+#qed_bw.tab_dndt_chi_max = 1000.0
+#qed_bw.tab_dndt_how_many = 256
+#qed_bw.tab_pair_chi_min = 0.01
+#qed_bw.tab_pair_chi_max = 1000.0
+#qed_bw.tab_pair_chi_how_many = 256
+#qed_bw.tab_pair_frac_how_many = 256
+#qed_bw.save_table_in = "bw_table"
+
+#qed_bw.lookup_table_mode = "load"
+#qed_bw.load_table_from = "bw_table"
+
+qed_qs.chi_min = 0.001
+
+qed_qs.lookup_table_mode = "builtin"
+
+qed_qs.photon_creation_energy_threshold = 0.0
+
+#qed_qs.lookup_table_mode = "generate"
+#qed_qs.tab_dndt_chi_min = 0.001
+#qed_qs.tab_dndt_chi_max = 1000.0
+#qed_qs.tab_dndt_how_many = 256
+#qed_qs.tab_em_chi_min = 0.001
+#qed_qs.tab_em_frac_min = 1.0e-12
+#qed_qs.tab_em_chi_max = 1000.0
+#qed_qs.tab_em_chi_how_many = 256
+#qed_qs.tab_em_frac_how_many = 256
+#qed_qs.save_table_in = "qs_table"
+
+#qed_qs.lookup_table_mode = "load"
+#qed_qs.load_table_from = "qs_table"
+#################################
+
+### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] )
+particles.E_ext_particle_init_style = "constant"
+particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863
+####
+
+### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] )
+particles.B_ext_particle_init_style = "constant"
+particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857
+####
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.period = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex
+diag1.p1.variables = ux uy uz w
+diag1.p2.variables = ux uy uz w
+diag1.p3.variables = ux uy uz w
+diag1.p4.variables = ux uy uz w
+
+diag1.qsp_1.variables = ux uy uz w
+diag1.qsp_2.variables = ux uy uz w
+diag1.qsp_3.variables = ux uy uz w
+diag1.qsp_4.variables = ux uy uz w
diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs_2d_tau_init b/Examples/Modules/qed/quantum_synchrotron/inputs_2d_tau_init
deleted file mode 100644
index a29adf287..000000000
--- a/Examples/Modules/qed/quantum_synchrotron/inputs_2d_tau_init
+++ /dev/null
@@ -1,122 +0,0 @@
-#################################
-####### GENERAL PARAMETERS ######
-#################################
-max_step = 1
-amr.n_cell =  128 128
-amr.max_grid_size = 128   # maximum size of each AMReX box, used to decompose the domain
-amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain
-geometry.coord_sys   = 0                  # 0: Cartesian
-geometry.is_periodic = 0     0      # Is periodic?
-geometry.prob_lo     =  -32.e-6   -32.e-6    # physical domain
-geometry.prob_hi     =  32.e-6   32.e-6
-amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
-
-#################################
-############ NUMERICS ###########
-#################################
-algo.current_deposition = esirkepov
-algo.charge_deposition = standard
-algo.field_gathering = energy-conserving
-algo.particle_pusher = boris
-interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
-interpolation.noy = 3
-interpolation.noz = 3
-warpx.verbose = 1
-warpx.do_dive_cleaning = 0
-warpx.use_filter = 1
-warpx.cfl = 1. # if 1., the time step is set to its CFL limit
-warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
-warpx.serialize_ics = 1
-
-
-#################################
-############ PLASMA #############
-#################################
-particles.species_names = electrons positrons qs_phot
-particles.photon_species = qs_phot
-#################################
-
-electrons.species_type = "electron"
-electrons.injection_style = "NUniformPerCell"
-electrons.profile = "constant"
-electrons.xmin = -30e-6
-electrons.ymin = -30e-6
-electrons.zmin = -30e-6
-electrons.xmax = 30e-6
-electrons.ymax = 30e-6
-electrons.zmax = 30e-6
-electrons.num_particles_per_cell_each_dim = 2 2
-electrons.density = 1e19
-electrons.profile = "constant"
-electrons.momentum_distribution_type = "gaussian"
-electrons.ux_m = 0.0
-electrons.uy_m = 0.0
-electrons.uz_m = 0.0
-electrons.ux_th = 100.
-electrons.uy_th = 100.
-electrons.uz_th = 100.
-##########QED####################
-electrons.do_qed = 1
-electrons.do_qed_quantum_sync = 1
-electrons.qed_quantum_sync_phot_product_species = qs_phot
-#################################
-
-positrons.species_type = "positron"
-positrons.injection_style = "NUniformPerCell"
-positrons.profile = "constant"
-positrons.xmin = -30e-6
-positrons.ymin = -30e-6
-positrons.zmin = -30e-6
-positrons.xmax = 30e-6
-positrons.ymax = 30e-6
-positrons.zmax = 30e-6
-positrons.num_particles_per_cell_each_dim = 2 2
-positrons.density = 1e19
-positrons.profile = "constant"
-positrons.momentum_distribution_type = "gaussian"
-positrons.ux_m = 0.0
-positrons.uy_m = 0.0
-positrons.uz_m = 0.0
-positrons.ux_th = 100.
-positrons.uy_th = 100.
-positrons.uz_th = 100.
-##########QED####################
-positrons.do_qed = 1
-positrons.do_qed_quantum_sync = 1
-positrons.qed_quantum_sync_phot_product_species = qs_phot
-#################################
-
-### PRODUCT SPECIES ###
-qs_phot.species_type = "photon"
-qs_phot.injection_style = nuniformpercell
-qs_phot.num_particles_per_cell_each_dim = 1 1
-qs_phot.profile = constant
-qs_phot.xmin = 1 ## Ugly trick to avoid photons at T=0
-qs_phot.xmax = -1 ##  Ugly trick to avoid photons at T=0
-qs_phot.density = 0.0
-qs_phot.momentum_distribution_type = "gaussian"
-#################################
-
-#######QED ENGINE PARAMETERS#####
-qed_qs.lookup_table_mode = "dummy_builtin"
-qed_qs.photon_creation_energy_threshold = 1.63742112993e-13
-
-#qed_qs.lookup_table_mode = "generate"
-#qed_qs.chi_min = 0.001
-#qed_qs.tab_dndt_chi_min = 0.001
-#qed_qs.tab_dndt_chi_max = 200
-#qed_qs.tab_dndt_how_many = 32
-#qed_qs.tab_em_chi_min = 0.001
-#qed_qs.tab_em_chi_max = 200
-#qed_qs.tab_em_chi_how_many = 2
-#qed_qs.tab_em_prob_how_many = 2
-#qed_qs.save_table_in = "qs_micro_table"
-
-#qed_qs.lookup_table_mode = "load"
-#qed_qs.load_table_from = "qs_micro_table"
-#################################
-
-# Diagnostics
-diagnostics.diags_names = diag1
-diag1.period = 1
-diag1.diag_type = Full
diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs_3d b/Examples/Modules/qed/quantum_synchrotron/inputs_3d
new file mode 100644
index 000000000..66670b4db
--- /dev/null
+++ b/Examples/Modules/qed/quantum_synchrotron/inputs_3d
@@ -0,0 +1,224 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell =  16 16 16
+amr.max_grid_size = 16   # maximum size of each AMReX box, used to decompose the domain
+amr.blocking_factor = 8 # minimum size of each AMReX box, used to decompose the domain
+geometry.coord_sys   = 0                  # 0: Cartesian
+geometry.is_periodic = 1  1  1      # Is periodic?
+geometry.prob_lo     =  -0.25e-6   -0.25e-6  -0.25e-6  # physical domain
+geometry.prob_hi     =  0.25e-6   0.25e-6   0.25e-6
+amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported)
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.current_deposition = esirkepov
+algo.charge_deposition = standard
+algo.field_gathering = energy-conserving
+algo.particle_pusher = boris
+interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z
+interpolation.noy = 3
+interpolation.noz = 3
+warpx.verbose = 1
+warpx.do_dive_cleaning = 0
+warpx.use_filter = 1
+warpx.cfl = 1. # if 1., the time step is set to its CFL limit
+warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition
+warpx.serialize_ics = 1
+
+#################################
+############ PLASMA #############
+#################################
+particles.species_names = p1 p2 p3 p4 qsp_1 qsp_2 qsp_3 qsp_4 dummy_ele dummy_pos
+particles.photon_species = qsp_1 qsp_2 qsp_3 qsp_4
+#################################
+
+p1.species_type = "electron"
+p1.injection_style = "NUniformPerCell"
+p1.profile = "constant"
+p1.num_particles_per_cell_each_dim = 8 8 4
+p1.density = 1e1
+p1.profile = "constant"
+p1.momentum_distribution_type = "gaussian"
+p1.ux_m = 10.0
+##########QED####################
+p1.do_qed = 1
+p1.do_qed_quantum_sync = 1
+p1.qed_quantum_sync_phot_product_species = qsp_1
+#################################
+
+p2.species_type = "electron"
+p2.injection_style = "NUniformPerCell"
+p2.profile = "constant"
+p2.num_particles_per_cell_each_dim = 8 8 4
+p2.density = 1e1
+p2.profile = "constant"
+p2.momentum_distribution_type = "gaussian"
+p2.uy_m = 100.0
+##########QED####################
+p2.do_qed = 1
+p2.do_qed_quantum_sync = 1
+p2.qed_quantum_sync_phot_product_species = qsp_2
+#################################
+
+p3.species_type = "positron"
+p3.injection_style = "NUniformPerCell"
+p3.profile = "constant"
+p3.num_particles_per_cell_each_dim = 8 8 4
+p3.density = 1e1
+p3.profile = "constant"
+p3.momentum_distribution_type = "gaussian"
+p3.uz_m = 1000.0
+##########QED####################
+p3.do_qed = 1
+p3.do_qed_quantum_sync = 1
+p3.qed_quantum_sync_phot_product_species = qsp_3
+#################################
+
+p4.species_type = "positron"
+p4.injection_style = "NUniformPerCell"
+p4.profile = "constant"
+p4.num_particles_per_cell_each_dim = 8 8 4
+p4.density = 1e1
+p4.profile = "constant"
+p4.momentum_distribution_type = "gaussian"
+p4.ux_m = 5773.502691896
+p4.uy_m = 5773.502691896
+p4.uz_m = 5773.502691896
+##########QED####################
+p4.do_qed = 1
+p4.do_qed_quantum_sync = 1
+p4.qed_quantum_sync_phot_product_species = qsp_4
+#################################
+
+### PRODUCT SPECIES ###
+qsp_1.species_type = "photon"
+qsp_1.injection_style = nuniformpercell
+qsp_1.num_particles_per_cell_each_dim = 0 0
+qsp_1.profile = constant
+qsp_1.density = 0.0
+qsp_1.momentum_distribution_type = "gaussian"
+qsp_1.do_qed = 1
+qsp_1.do_qed_breit_wheeler = 1
+qsp_1.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_1.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_2.species_type = "photon"
+qsp_2.injection_style = nuniformpercell
+qsp_2.num_particles_per_cell_each_dim = 0 0
+qsp_2.profile = constant
+qsp_2.density = 0.0
+qsp_2.momentum_distribution_type = "gaussian"
+qsp_2.do_qed = 1
+qsp_2.do_qed_breit_wheeler = 1
+qsp_2.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_2.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_3.species_type = "photon"
+qsp_3.injection_style = nuniformpercell
+qsp_3.num_particles_per_cell_each_dim = 0 0
+qsp_3.profile = constant
+qsp_3.density = 0.0
+qsp_3.momentum_distribution_type = "gaussian"
+qsp_3.do_qed = 1
+qsp_3.do_qed_breit_wheeler = 1
+qsp_3.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_3.qed_breit_wheeler_pos_product_species = dummy_pos
+
+qsp_4.species_type = "photon"
+qsp_4.injection_style = nuniformpercell
+qsp_4.num_particles_per_cell_each_dim = 0 0
+qsp_4.profile = constant
+qsp_4.density = 0.0
+qsp_4.momentum_distribution_type = "gaussian"
+qsp_4.do_qed = 1
+qsp_4.do_qed_breit_wheeler = 1
+qsp_4.qed_breit_wheeler_ele_product_species = dummy_ele
+qsp_4.qed_breit_wheeler_pos_product_species = dummy_pos
+
+#################################
+
+dummy_ele.species_type = "electron"
+dummy_ele.injection_style = nuniformpercell
+dummy_ele.num_particles_per_cell_each_dim = 0 0
+dummy_ele.profile = constant
+dummy_ele.density = 0.0
+dummy_ele.momentum_distribution_type = "gaussian"
+dummy_ele.do_qed = 0
+
+dummy_pos.species_type = "positron"
+dummy_pos.injection_style = nuniformpercell
+dummy_pos.num_particles_per_cell_each_dim = 0 0
+dummy_pos.profile = constant
+dummy_pos.density = 0.0
+dummy_pos.momentum_distribution_type = "gaussian"
+dummy_pos.do_qed = 0
+
+
+#################################
+
+##########QED TABLES####################
+qed_bw.chi_min = 0.001
+
+qed_bw.lookup_table_mode = "builtin"
+
+#qed_bw.lookup_table_mode = "generate"
+#qed_bw.tab_dndt_chi_min = 0.01
+#qed_bw.tab_dndt_chi_max = 1000.0
+#qed_bw.tab_dndt_how_many = 256
+#qed_bw.tab_pair_chi_min = 0.01
+#qed_bw.tab_pair_chi_max = 1000.0
+#qed_bw.tab_pair_chi_how_many = 256
+#qed_bw.tab_pair_frac_how_many = 256
+#qed_bw.save_table_in = "bw_table"
+
+#qed_bw.lookup_table_mode = "load"
+#qed_bw.load_table_from = "bw_table"
+
+qed_qs.chi_min = 0.001
+
+qed_qs.lookup_table_mode = "builtin"
+
+qed_qs.photon_creation_energy_threshold = 0.0
+
+#qed_qs.lookup_table_mode = "generate"
+#qed_qs.tab_dndt_chi_min = 0.001
+#qed_qs.tab_dndt_chi_max = 1000.0
+#qed_qs.tab_dndt_how_many = 256
+#qed_qs.tab_em_chi_min = 0.001
+#qed_qs.tab_em_frac_min = 1.0e-12
+#qed_qs.tab_em_chi_max = 1000.0
+#qed_qs.tab_em_chi_how_many = 256
+#qed_qs.tab_em_frac_how_many = 256
+#qed_qs.save_table_in = "qs_table"
+
+#qed_qs.lookup_table_mode = "load"
+#qed_qs.load_table_from = "qs_table"
+#################################
+
+### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] )
+particles.E_ext_particle_init_style = "constant"
+particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863
+####
+
+### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] )
+particles.B_ext_particle_init_style = "constant"
+particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857
+####
+
+# Diagnostics
+diagnostics.diags_names = diag1
+diag1.period = 1
+diag1.diag_type = Full
+diag1.fields_to_plot = Ex
+diag1.p1.variables = ux uy uz w
+diag1.p2.variables = ux uy uz w
+diag1.p3.variables = ux uy uz w
+diag1.p4.variables = ux uy uz w
+
+diag1.qsp_1.variables = ux uy uz w
+diag1.qsp_2.variables = ux uy uz w
+diag1.qsp_3.variables = ux uy uz w
+diag1.qsp_4.variables = ux uy uz w
diff --git a/Examples/Modules/qed/schwinger/analysis_schwinger.py b/Examples/Modules/qed/schwinger/analysis_schwinger.py
index f52512478..25dadeb42 100755
--- a/Examples/Modules/qed/schwinger/analysis_schwinger.py
+++ b/Examples/Modules/qed/schwinger/analysis_schwinger.py
@@ -21,9 +21,9 @@ import checksumAPI
 # define some parameters
 
 c = 299792458.
-m_e = 9.10938356e-31
-e = 1.6021766208e-19
-hbar = 1.054571800e-34
+m_e = 9.1093837015e-31
+e =1.602176634e-19
+hbar = 1.054571817e-34
 E_S = m_e**2*c**3/e/hbar # Schwinger field
 
 dV = (1.e-6)**3 # total simulation volume