diff --git a/docs/changes/275.maintenance.rst b/docs/changes/275.maintenance.rst
new file mode 100644
index 000000000..71ce09633
--- /dev/null
+++ b/docs/changes/275.maintenance.rst
@@ -0,0 +1 @@
+Added filling of CREF keyword (IRF axis order) in the output files
diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py
index dd0d22f32..e14c20670 100644
--- a/pyirf/io/gadf.py
+++ b/pyirf/io/gadf.py
@@ -76,6 +76,8 @@ def create_aeff2d_hdu(
     header["HDUCLAS3"] = "POINT-LIKE" if point_like else "FULL-ENCLOSURE"
     header["HDUCLAS4"] = "AEFF_2D"
     header["DATE"] = Time.now().utc.iso
+    idx = aeff.colnames.index("EFFAREA") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,THETA_LO:THETA_HI)"
     _add_header_cards(header, **header_cards)
 
     return BinTableHDU(aeff, header=header, name=extname)
@@ -133,6 +135,8 @@ def create_psf_table_hdu(
     header["HDUCLAS3"] = "FULL-ENCLOSURE"
     header["HDUCLAS4"] = "PSF_TABLE"
     header["DATE"] = Time.now().utc.iso
+    idx = psf_.colnames.index("RPSF") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,THETA_LO:THETA_HI,RAD_LO:RAD_HI)"
     _add_header_cards(header, **header_cards)
 
     return BinTableHDU(psf_, header=header, name=extname)
@@ -191,6 +195,8 @@ def create_energy_dispersion_hdu(
     header["HDUCLAS3"] = "POINT-LIKE" if point_like else "FULL-ENCLOSURE"
     header["HDUCLAS4"] = "EDISP_2D"
     header["DATE"] = Time.now().utc.iso
+    idx = edisp.colnames.index("MATRIX") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,MIGRA_LO:MIGRA_HI,THETA_LO:THETA_HI)"
     _add_header_cards(header, **header_cards)
 
     return BinTableHDU(edisp, header=header, name=extname)
@@ -247,6 +253,8 @@ def create_background_2d_hdu(
     header["HDUCLAS3"] = "FULL-ENCLOSURE"
     header["HDUCLAS4"] = "BKG_2D"
     header["DATE"] = Time.now().utc.iso
+    idx = bkg.colnames.index("BKG") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,THETA_LO:THETA_HI)"
     _add_header_cards(header, **header_cards)
 
     return BinTableHDU(bkg, header=header, name=extname)
@@ -300,6 +308,8 @@ def create_rad_max_hdu(
     header["HDUCLAS3"] = "POINT-LIKE"
     header["HDUCLAS4"] = "RAD_MAX_2D"
     header["DATE"] = Time.now().utc.iso
+    idx = rad_max_table.colnames.index("RAD_MAX") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,THETA_LO:THETA_HI)"
     _add_header_cards(header, **header_cards)
 
     return BinTableHDU(rad_max_table, header=header, name=extname)
diff --git a/pyirf/io/tests/test_gadf.py b/pyirf/io/tests/test_gadf.py
index 62139af3a..653834d8b 100644
--- a/pyirf/io/tests/test_gadf.py
+++ b/pyirf/io/tests/test_gadf.py
@@ -181,3 +181,32 @@ def test_rad_max_schema(rad_max_hdu):
 
     _, hdu = rad_max_hdu
     RAD_MAX.validate_hdu(hdu)
+
+
+def test_cref(aeff2d_hdus, edisp_hdus, psf_hdu, bg_hdu, rad_max_hdu):
+    for point_like in range (2):
+        hdus = [fits.PrimaryHDU(),
+                aeff2d_hdus[1][point_like],
+                edisp_hdus[1][point_like],
+                psf_hdu[1],
+                bg_hdu[1],
+                rad_max_hdu[1]]
+
+        with tempfile.NamedTemporaryFile(suffix='.fits') as f:
+            fits.HDUList(hdus).writeto(f.name)
+            readhdul = fits.open(f.name)
+            for hdu in readhdul[1:]: # skip Primary
+                names=hdu.data.columns.names
+                # the IRF is always in the last column
+                cref = hdu.header['CREF'+str(len(names))]
+                # drop parentheses and get the expected list of axes
+                cref = cref[1:-1].split(',')
+                # transposed due to different fits and numpy axis order
+                readshape=hdu.data[names[-1]].T.shape
+                # only one set of IRFs saved per test file
+                assert readshape[-1] == 1, "first axis size is not 1"
+
+                for i in range(len(readshape)-1):
+                    for axis in cref[i].split(':'):
+                        err="not matching shape in "+hdu.header['EXTNAME']+" at axis "+axis
+                        assert readshape[i] == hdu.data[axis].shape[1], err