Crab Nebula Performance

Crab Nebula Performance

Draft versions

Version 1, 20220923: [[1]]

Version 2, 20230124: [[2]]

Version 3, 20230126: [[3]]

Version 4, 20230321: [[4]]

Version 5, 20230329: [[5]]

Format was changes after A&A considered the paper topic was not suitable for their journal

Version 6, sent to ApJ, 20230421: [[6]]

Version 9, addressing the ApJ referee's comments, 20230623: [[7]]

Version 10, minor changes (acknowledgements, author list, fig 16 typo, zenodo doi for repo containing the paper's plots data), 20230719: [[8]]

General information

• Name of the source: Crab Nebula and pulsar
• Brief description of the source:
  • Object type: PWN
  • Other relevant information
  • RA: 05 34 31.9 (hh mm ss), Dec: +20 00 52.2 (dd mm ss)
  • RA, Dec in deg (ICRS): 83.6331, 22.0145

• Analysis by:
  • Abelardo Moralejo (IFAE, moralejo@ifae.es)
  • Rubén López-Coto (IAA-CSIC, rlopezcoto@iaa.es)
  • Thomas Vuillaume (LAPP, thomas.vuillaume@lapp.in2p3.fr)
  • Seiya Nozaki (MPP, nozaki@mpp.mpg.de)
  • Daniel Morcuende (UCM, now at IAA-CSIC, dmorcuende@iaa.es)

Data-taking information

• Data selection based on the notebook: https://github.com/cta-observatory/cta-lstchain/blob/master/notebooks/data_selection.ipynb
  • Excluded runs with issues in interleaved pedestals
  • Run selection also based on "healthy intensity spectra", i.e., rate of cosmics with intensity between 80 and 120 p.e. larger than 800 events/s.
  • Runs taken in wobble mode. Selected wobble pointings at ~0.4 deg from the source [0.35, 0.45] deg

• Observation time: 35.94 hours

• Dates of data-taking:

1 : 2020-11-17 : [2914]
2 : 2020-11-18 : [2929, 2930, 2931, 2932, 2933, 2934]
3 : 2020-11-19 : [2949, 2950]
4 : 2020-11-20 : [2967, 2968, 2969, 2970, 2971, 2972, 2973, 2974, 2975, 2976, 2977]
5 : 2020-11-21 : [2988, 2989, 2990, 2991, 2992]
6 : 2020-11-22 : [3004, 3005, 3006, 3007, 3008]
7 : 2020-12-07 : [3093, 3094, 3095, 3096]
8 : 2020-12-13 : [3231, 3232]
9 : 2020-12-14 : [3243]
10 : 2020-12-15 : [3270, 3271, 3272, 3273, 3274, 3275, 3276, 3277, 3278, 3279]
11 : 2020-12-17 : [3318, 3319, 3320, 3321]
12 : 2020-12-18 : [3328, 3329, 3330]
13 : 2020-12-19 : [3338, 3339, 3340]
14 : 2020-12-20 : [3355, 3356]
15 : 2020-12-21 : [3373]
16 : 2021-02-11 : [3598, 3599, 3600, 3601]
17 : 2021-02-12 : [3615]
18 : 2021-02-13 : [3632, 3633, 3634, 3635]
19 : 2021-02-15 : [3672, 3673, 3674, 3675, 3676, 3677]
20 : 2021-02-16 : [3706, 3707, 3708]
21 : 2021-03-15 : [4067, 4068]
22 : 2021-03-16 : [4086, 4087]
23 : 2021-09-04 : [6045]
24 : 2021-09-05 : [6073]
25 : 2021-09-14 : [6304]
26 : 2022-02-04 : [6872, 6873, 6874, 6875]
27 : 2022-02-05 : [6892, 6893, 6894, 6895]
28 : 2022-02-23 : [7097, 7098, 7099]
29 : 2022-02-28 : [7133, 7136]
30 : 2022-03-01 : [7161]
31 : 2022-03-02 : [7195, 7196, 7197, 7199, 7200]
32 : 2022-03-03 : [7227, 7228, 7231, 7232, 7233]
33 : 2022-03-04 : [7253, 7254, 7255, 7256]
34 : 2022-03-05 : [7274, 7275, 7276, 7277]

• List of all selected runs:

[2914, 2929, 2930, 2931, 2932, 2933, 2934, 2949, 2950, 2967, 2968, 2969,
 2970, 2971, 2972, 2973, 2974, 2975, 2976, 2977, 2988, 2989, 2990, 2991,
 2992, 3004, 3005, 3006, 3007, 3008, 3093, 3094, 3095, 3096, 3231, 3232,
 3243, 3270, 3271, 3272, 3273, 3274, 3275, 3276, 3277, 3278, 3279, 3318,
 3319, 3320, 3321, 3328, 3329, 3330, 3338, 3339, 3340, 3355, 3356, 3373,
 3598, 3599, 3600, 3601, 3615, 3632, 3633, 3634, 3635, 3672, 3673, 3674,
 3675, 3676, 3677, 3706, 3707, 3708, 4067, 4068, 4086, 4087, 6045, 6073,
 6304, 6872, 6873, 6874, 6875, 6892, 6893, 6894, 6895, 7097, 7098, 7099,
 7133, 7136, 7161, 7195, 7196, 7197, 7199, 7200, 7227, 7228, 7231, 7232,
 7233, 7253, 7254, 7255, 7256, 7274, 7275, 7276, 7277]

• ZD range (deg): 0 - 35 deg
• AZ range (deg): 90 - 270 deg

• Transmission at 9 km: See corresponding values in the ELOG entries. We do not use it for data selection. Instead we apply selection cuts based on rates (see above).
• Joint observations with MAGIC? : No

Monte Carlo information

• Link to MC files used:
  • Particle types:
  • ZD (deg): Zenith Range
  • AZ (deg): Azimuth Range

Pointing directions of the nodes where MC training and testing samples were simualted.

Analysis Code

The code used in the analysis can be found at https://github.com/vuillaut/crab-performance-paper or in the ZIP archive there:

/fefs/aswg/software/lst-crab-performance-paper/crab-performance-paper-1.0.zip

DL1 data

Used standard DL1b files produced by lstosa (using lstchain v0.9.x):

• real data: tailcut8-4 (plus cleaning based on pixel-wise pedestal std dev & dynamic cleaning, see paper):

/fefs/aswg/data/real/DL1/YYYYMMDD/v0.9/tailcut84/dl1_LST-1.Run?????.h5

• MC: tailcut8-4 (plus cleaning based on pixel-wise pedestal std dev, dynamic cleaning). NSB tuning to Crab:

/fefs/aswg/data/mc/DL1/AllSky/20221027_v0.9.9_crab_tuned/TrainingDataset/dec_2276/
/fefs/aswg/data/mc/DL1/AllSky/20221027_v0.9.9_crab_tuned/TestingDataset/ (see paper, closest pointing node used for each data run)

Random forest

Please include your specific .json files used in producing your RFs and for which source analysis (source dependent or independent), and variabules used in the RF.

• lstchain-0.9.9
• source-independent:

config: /fefs/aswg/workspace/seiya.nozaki/data/MC/v0.7.5/tailcut84_dynamic_bllac/srcdep/RF/lstchain_src_dep_config.json
models: /fefs/aswg/data/models/AllSky/20221027_v0.9.9_crab_tuned/dec_2276/
lstmcpipe configuration: https://github.com/cta-observatory/lstmcpipe/tree/master/production_configs/20221027_v0.9.9_crab_tuned

• source-dependent:

config: /fefs/aswg/workspace/seiya.nozaki/Crab_performance_paper/20221027_v0.9.9_crab_tuned/combined_off_axis_1deg/RF/lstchain_src_dep_config.json
models: /fefs/aswg/workspace/seiya.nozaki/Crab_performance_paper/20221027_v0.9.9_crab_tuned/combined_off_axis_1deg/RF

DL2 data

Information about your DL2 data and settings such as: specific versions of lstchain, specific .json version used.

• lstchain-0.9.9
• source-indep

/fefs/aswg/workspace/abelardo.moralejo/Crab_performance_paper/data_v0.9.9/DL2/dl2*.h5  (selected runs for final analysis)

Closest MC node to each data run:

• theta_10.0_az_102.199:

2929, 2930, 2931, 2950, 2970, 2971, 2972, 2973, 2988,
3095, 3096, 3273, 3274, 3275, 3318, 3319, 3320, 3338

• theta_10.0_az_248.117

2914, 2932, 2933, 2934, 2974, 2975, 2976, 2989, 2990,                                                                                              
2991, 3004, 3005, 3006, 3276, 3277, 3278, 3279, 3321,                                                                                              
3328, 3329, 3330, 3339, 3340, 3355, 3356, 3598, 3632,                                                                                              
3672, 3673, 3674, 3675, 6872, 6873, 6874, 6892, 6893,                                                                                              
7133, 7161, 7195, 7196, 7197, 7227, 7228, 7253, 7274

• theta_23.63_az_100.758

2949, 2968, 2969, 3093, 3094, 3271, 3272, 6304

• theta_23.63_az_259.265

2977, 2992, 3007, 3008, 3231, 3243, 3373, 3599,                                                                                                    
3600, 3615, 3633, 3634, 3676, 3677, 3706, 3707,                                                                                                    
4067, 4086, 6875, 6894, 7097, 7098, 7099, 7199,                                                                                                    
7231, 7254, 7255, 7256, 7275

• theta_32.059_az_102.217

2967, 3270, 6045, 6073

• theta_32.059_az_248.099

3232, 3601, 3635, 3708, 4068, 4087, 6895, 7136,                                                                                                    
7200, 7232, 7233, 7276, 7277 

Closest MC pointing node to each Crab data run.

• source-dep

real data: /fefs/aswg/workspace/seiya.nozaki/Crab_performance_paper/20221027_v0.9.9_crab_tuned/combined_off_axis_1deg/DL2/data/dl2*.h5
MC data: /fefs/aswg/workspace/seiya.nozaki/Crab_performance_paper/20221027_v0.9.9_crab_tuned/combined_off_axis_1deg/DL2_MC/data

IRF creation

For the production of IRFs for high-level analysis, the ones that are included in the DL3 files, we used the script lstchain_create_irf_files.

Input test DL2 MC used:

/fefs/aswg/data/mc/DL2/AllSky/20221027_v0.9.9_crab_tuned/TestingDataset/dec_2276

IRF files were created for each MC pointing node used for the analysis of the Crab data runs corresponding to the following subdirectories:

theta_10.0_az_102.199
theta_10.0_az_248.117
theta_23.63_az_100.758
theta_23.63_az_259.265
theta_32.059_az_102.217
theta_32.059_az_248.099

Command line options used:

point-like                                                                                                                                                                                                                                                                                                                        
energy-dependent-gh                                                                                                                                                            
energy-dependent-theta

For gh-efficiency and theta-containment, we used a grid of values: [0.4, 0.5, 0.6, 0.7, 0.8, 0.9] for gh-efficiency (efficiency of gammaness cut) and [0.7, 0.8, 0.9] for theta-containment (efficiency of theta cut), in all possible combinations.

Config file used (the same used for the creation of the DL3 files described in the next section):

{
  "EventSelector": {
    "filters": {
      "intensity": [80, Infinity],
      "width": [0, Infinity],
      "length": [0, Infinity],
      "r": [0, 1],
      "wl": [0, 1],
      "leakage_intensity_width_2": [0, 1],
      "event_type": [32, 32]
    }
  },
  "DL3Cuts": {   
    "min_event_p_en_bin": 100,
    "min_gh_cut": 0.1,
    "max_gh_cut": 0.95,
    "min_theta_cut": 0.05,
    "max_theta_cut": 0.32,
    "fill_theta_cut": 0.32,
    "allowed_tels": [1]
  },
  "DataBinning": {
    "true_energy_min": 0.002,
    "true_energy_max": 200,
    "true_energy_n_bins_per_decade": 5,
    "reco_energy_min": 0.002,
    "reco_energy_max": 200,
    "reco_energy_n_bins_per_decade": 5,
    "energy_migration_min": 0.2,
    "energy_migration_max": 5,
    "energy_migration_n_bins": 31,
    "fov_offset_min": 0.1,
    "fov_offset_max": 1.1,
    "fov_offset_n_edges": 9,
    "bkg_fov_offset_min": 0,
    "bkg_fov_offset_max": 10,
    "bkg_fov_offset_n_edges": 21,
    "source_offset_min": 0.0001,
    "source_offset_max": 1.0001,
    "source_offset_n_edges": 1000
  }
}

DL3 data selection

Information about DL3 data selection.

• source-independent analysis

- For the analysis of the entire dataset from Nov 2020 to Mar 2022 we applied an intensity cut of >80 p.e. - For the analysis of the dataset >Aug 2021 (when trigger settings became stable) we applied an intensity cut of >50 p.e.

We did not apply further cuts in any other parameter (other than intensity, gammaness and theta).

The config file used to produce the DL3 and corresponding IRFs is:

{
  "EventSelector": {
    "filters": {
      "intensity": [80, Infinity],
      "width": [0, Infinity],
      "length": [0, Infinity],
      "r": [0, 1],
      "wl": [0, 1],
      "leakage_intensity_width_2": [0, 1],
      "event_type": [32, 32]
    }
  },
  "DL3Cuts": {   
    "min_event_p_en_bin": 100,
    "min_gh_cut": 0.1,
    "max_gh_cut": 0.95,
    "min_theta_cut": 0.05,
    "max_theta_cut": 0.32,
    "fill_theta_cut": 0.32,
    "allowed_tels": [1]
  },
  "DataBinning": {
    "true_energy_min": 0.002,
    "true_energy_max": 200,
    "true_energy_n_bins_per_decade": 5,
    "reco_energy_min": 0.002,
    "reco_energy_max": 200,
    "reco_energy_n_bins_per_decade": 5,
    "energy_migration_min": 0.2,
    "energy_migration_max": 5,
    "energy_migration_n_bins": 31,
    "fov_offset_min": 0.1,
    "fov_offset_max": 1.1,
    "fov_offset_n_edges": 9,
    "bkg_fov_offset_min": 0,
    "bkg_fov_offset_max": 10,
    "bkg_fov_offset_n_edges": 21,
    "source_offset_min": 0.0001,
    "source_offset_max": 1.0001,
    "source_offset_n_edges": 1000
  }
}

After the pre-selection of events based on intensity, the final selection of gamma-ray-like events is based on energy-dependent gammaness and theta selection cuts with a given efficiency. The baseline efficiency in the gamma-ray selection cuts is 70% for these two parameters. We also produce DL3 files assuming a grid of efficiency values: 40-90% for gammaness and 70-90% for theta. All these efficiency cut combinations are used to test the stability of the resulting spectrum against the choice of the cut values (see Spectral results).

• source-dependent analysis

min_event_p_en_bin": 100,
min_gh_cut": 0.1,
max_gh_cut": 0.95,
min_alpha_cut": 1,
max_alpha_cut": 20,
fill_alpha_cut": 20,

DL3 data

• source-independent analysis

/fefs/aswg/workspace/daniel.morcuende/data/real/DL3/Crab_performance/AllSkyMC_v0.9.9/intensity{50,80}/all_nodes/gh_eff_0.7_th_cont_0.7

In particular, the DL3 dataset used for calculating the baseline SED and light curve is:

/fefs/aswg/workspace/daniel.morcuende/data/real/DL3/Crab_performance/AllSkyMC_v0.9.9/intensity80/all_nodes/gh_eff_0.7_th_cont_0.7

The rest of the DL3 datasets with different efficiencies are used to test the stability of the SED against the change in the efficiency of selection cuts.

Warning: These DL3 files were produced with a buggy version of the pyIRF code with a wrong normalization of EDISP (see issue #249 in pyirf)

High-level analysis

Information about the production of higher-level analysis.

• lstchain 0.9.9 to generate source-indep, source-dep IRFs and DL3 files (see DL3 data)
  • Baseline analysis uses the DL3 sample obtained assuming 70% efficiency in gammaness and theta gamma-ray selection cuts
  • Note: It is recommended to produce the DL3 files unzipped. Otherwise, loading the files with Gammapy could take a substantially larger run time in the unzipping part (10 times more than with the unzipped DL3 files).
• Science Tool: Gammapy 0.20.1
• Point-like IRF (closest MC pointing node), 1D ON-OFF analysis
• Livetime (effective time) of the entire observation dataset (Nov 2020 - Mar 2022): 34.2 h

Gammapy settings

Source-independent analysis

As the DL3 files were produced with energy-dependent theta cuts, we need to include RAD_MAX IRF field when inputting the observations:

required_irf=["aeff", "edisp", "rad_max"]

• ON region is defined using the PointSkyRegion class.

• We use a single OFF reflected region defined as:

region_finder = WobbleRegionsFinder(n_off_regions=1)
bkg_maker = ReflectedRegionsBackgroundMaker(region_finder=region_finder)

• True and estimated energy axes:

energy_axis = MapAxis.from_energy_bounds(
    0.01, 100, nbin=8, per_decade=True, unit="TeV", name="energy"
)
energy_axis_true = MapAxis.from_energy_bounds(
    0.005, 200, nbin=20, per_decade=True, unit="TeV", name="energy_true"
)

• We do not use any safe mask when creating the observation's Dataset object.

• We followed the stacked-analysis approach after data reduction from DL3 to DL4:

stacked_dataset = datasets.stack_reduce(name="LST-1")

• Energy range for the spectral fitting (50 GeV to 30 TeV):

energy_range = [0.05, 30] * u.TeV
stacked_dataset.mask_fit = stacked.counts.geom.energy_mask(energy_range[0], energy_range[1])

• For the joint analysis Fermi-LAT & LST-1, we only consider Fermi-LAT data points above 2 GeV and LST-1 SED up to 2 TeV.

• For the test of assuming +1% in the background normalization for the SED calculation we use:

stacked_dataset.acceptance_off *= 0.99

Analysis Results

Higher-level analysis results included in the publication.

Spectral results

• LST-1 SED model (log-parabolic parametrization) from 50 GeV to 30 TeV. It includes the statistical uncertainty band. Solid line and band.
• LST-1 & Fermi-LAT joint model (2 GeV to 2 TeV), also including the statistical uncertainty band. Dotted line.
• SED flux points (from a DL3 sample obtained assuming 70% efficiency in gammaness and theta gamma-ray selection cuts). Filled circles markers.
• Error band (hatched blue) encompassing statistical uncertainty bands from SEDs models based on datasets having different gamma-ray selection efficiencies. The efficiencies of the cuts were chosen in a grid of 40%-90% for gammaness and 70%-90% for angular cut (theta).
• Flux points testing a +1% change in the background normalization, which only affects the low-energy spectral points (light open markers).

For reference:

• Fermi-LAT Crab Nebula SED from Arakawa et al. (2020). Squared markers.
• MAGIC Crab Nebula SED from Aleksić et al. (2015). Dashed line.

Source-independent analysis

   Number of total bins            : 32 
   Number of fit bins              : 21 
   
   Fit statistic type              : wstat
   Fit statistic value (-2 log(L)) : 48.53
   
   Number of parameters            : 4
   Number of free parameters       : 3
   Chi2/ndf                        : 48.53/(21-3)
   
   Spectral model type       : LogParabolaSpectralModel
   Parameters:
     amplitude                     :   3.05e-10   +/- 2.2e-12 1 / (cm2 s TeV)
     reference             (frozen):      0.400       TeV         
     alpha                         :       2.25   +/-    0.01             
     beta                          :      0.114   +/-   0.006             
 

Comparison pre- and post-August 2021 (when trigger settings were frozen)

• Flux points from dataset after August 2021
  • 70% efficiency in gammaness and theta selection cuts
  • Intensity > 50 p.e.
• Best-fit model for dataset after August 2021
  • 70% efficiency in gammaness and theta selection cuts
  • Intensity > 50 p.e.
• Band from best-fit models taking into account several efficiencies in selection cuts for dataset after August 2021
  • Grid of selection cut efficiencies: 40-90% in gammaness and 70-90% in theta
  • Intensity > 50 p.e.
• Best-fit model for dataset before August 2021
  • 70% efficiency in gammaness and theta selection cuts
  • Intensity > 80 p.e.

Light curve results

* E > 100 GeV
* Assumed best-fit models calculated above
* Night-wise (34 nights)

Source-independent analysis

Weighted mean flux: 4.947e-10 +/- 3.059e-12 cm-2 s-1
Chi2: 119.2, ndf: 33, P-value: 1.11e-11
Compatible with constant flux hypothesis if assuming a systematic error on the nightly fluxes of 5.8% (added quadratically to statistical uncertainties) 
  -> Chi2: 32.7, ndf: 33, P-value: 4.80e-01 (i.e., compatible with constant flux hypothesis)

Pulsar analysis

For the Crab Pulsar analysis we applied the same cuts as those applied for source-dependent analysis, taken from:

/fefs/aswg/workspace/seiya.nozaki/Crab_performance_paper/20221027_v0.9.9_crab_tuned/combined_off_axis_1deg/DL3/8bin_80pe/g0.7_a0.7/irf/irf_20221027_v0.9.9_crab_tuned_node_theta_10.0_az_102.199__merged.fits.gz

The DL2 files with the phases added are in:

/fefs/aswg/workspace/ruben.lopez/analysis/crabnebula/data/2022_03_Performance_paper/dl2_src_dep_pulsar_parameters_20230206_test_Seiya_ring_wobble_up_to_1deg_extra_params/

The notebook to process them:

/fefs/aswg/workspace/ruben.lopez/analysis/crabnebula/data/2022_03_Performance_paper/notebooks/Crab_pulsar_src_dependent_efficiency_cuts_test_Seiya_ring_wobble_up_to_1deg_extra_params.ipynb

and to plot the results:

/fefs/aswg/workspace/ruben.lopez/analysis/crabnebula/data/2022_03_Performance_paper/notebooks/Plot_phaseogram.ipynb

The phases selected for P1, P2 and the OFF region to carry out the analysis are:

P1  = [0.983, 0.026]
P2  = [0.377, 0.422]
OFF = [0.520, 0.870]

Results after EDISP correction

Same analysis settings (v0.9), just correcting the EDISP normalization

LC

Average (weighted mean) flux >100 GeV:

5.348e-10 +/- 3.310e-12 cm-2 s-1
Chi2: 121.6, ndf: 33
P-value: 4.57e-12

SED

f0 = 3.26 +- 0.02 e-10 TeV-1 s-1 cm-2
alpha = 2.217 +- 0.009
beta = 0.137 +- 0.006
E0 = 400 GeV
Chi2/ndf = 58.4 / 18