Analysis April 2022 with lstchain v0.9.6 Pirola (cross-check)

From my_wiki
Revision as of 12:11, 11 July 2022 by Gpirola (talk | contribs) (2D Analysis)
Jump to: navigation, search

General information

  • Name of the source: LHAASO J2108+5157
  • Brief description of the source:
- Object type : Unidentified Galactic PeVatron candidate
- Distance (pc) : Unknown
- RA, Dec in deg (ICRS): 317.22, 51.95
  • Analysis by Giorgio Pirola (MPP, gpirola@mpp.mgp.de)

Data-taking information (Run selection)

  • Run selection with the use of the notebook from Abelardo Moralejo [1]
  • Atmospheric transmission extracted from ELOG
  • Zenith < 55 deg
  • Summary of selection cuts:
- Wobble in (0.45, 0.55)deg
- elapsed_time > 5 min
- transmission_cut > 0.65
- pedestal charge std dev < 1.8 p.e.
- Cosmic rate > 3000 ev/s
- Cosmic rate( > 10 p.e.) > 20 ev/s
- Cosmic rate( > 30 p.e.) > 3 ev/s
- muon ring width std dev < 0.023
  • Run selection after cuts:
 Total wobble runs: 177
 Observation time: 51.40 hours
 Selected Runs:
 1 : 2021-06-04 : [4913, 4914, 4915, 4916, 4917]
 2 : 2021-06-05 : [4935, 4936]
 3 : 2021-06-12 : [5028, 5029, 5030, 5031]
 4 : 2021-06-30 : [5071, 5072]
 5 : 2021-07-01 : [5080, 5081, 5082, 5083, 5084]
 6 : 2021-07-02 : [5091, 5092, 5093]
 7 : 2021-07-03 : [5101, 5102, 5103, 5104, 5105, 5106, 5107, 5108]
 8 : 2021-07-04 : [5115, 5116, 5117, 5118, 5119, 5120, 5121]
 9 : 2021-07-05 : [5135, 5136, 5137, 5138, 5139, 5140, 5141, 5142]
 10 : 2021-07-15 : [5270, 5272]
 11 : 2021-08-01 : [5440, 5441, 5442]
 12 : 2021-08-03 : [5461, 5462, 5463, 5464, 5465]
 13 : 2021-08-04 : [5473, 5474, 5475, 5476, 5477, 5478, 5479, 5480]
 14 : 2021-08-05 : [5491, 5492, 5493, 5494, 5497, 5498, 5499, 5500]
 15 : 2021-08-06 : [5505, 5506, 5507, 5508, 5509, 5510, 5511, 5512, 5513, 5514, 5515, 5516, 5517]
 16 : 2021-08-08 : [5576]
 17 : 2021-08-09 : [5590, 5591]
 18 : 2021-08-10 : [5641, 5642, 5643]
 19 : 2021-08-11 : [5681, 5682, 5683, 5684, 5685, 5686, 5687]
 20 : 2021-08-12 : [5707, 5708, 5709, 5710, 5711, 5712, 5713]
 21 : 2021-08-13 : [5727]
 22 : 2021-09-01 : [5947, 5948, 5949, 5950, 5952]
 23 : 2021-09-02 : [5980, 5981, 5982, 5983, 5984, 5985, 5986, 5987, 5988, 5989, 5990, 5991]
 24 : 2021-09-03 : [5999, 6000, 6001, 6002, 6003, 6004, 6005, 6006, 6007, 6008, 6009, 6010]
 25 : 2021-09-04 : [6023, 6024, 6034, 6035, 6036, 6037, 6038]
 26 : 2021-09-05 : [6057, 6058, 6059, 6060, 6061, 6062, 6063, 6064, 6065, 6066]
 27 : 2021-09-06 : [6079, 6080, 6082, 6083, 6084, 6085]
 28 : 2021-09-07 : [6130, 6131, 6132, 6133, 6134]
 29 : 2021-09-09 : [6175, 6176, 6177, 6178, 6179, 6180, 6181, 6182, 6183]
 30 : 2021-09-11 : [6230, 6231, 6233]
 31 : 2021-09-12 : [6254, 6255, 6256, 6257]

MC information

- AllSkyMC production
- standard DL1 produced in lstmcpipe v0.7.4 (Training DiffuseGammas and Protons, Testing Ring-like gammas), and v0.8.2 (Testing DiffuseGammas)
- RFs training data reconstruction up to DL2 done manually in lstchain v0.9.6 (lstmcpipe wasn't ready for the full AllSky MC at the time of analysis), IRFs created in dev lstchain (close to v0.9.7)
- For IRFs we used only the closest (in zenith) testing nodes to the LHAASO path


  • Nodes used for IRFs:
- node_theta_14.984_az_355.158_
- node_theta_32.059_az_355.158_
- node_theta_43.197_az_262.712_
- node_theta_52.374_az_301.217_

Lhaaso path delta nodes.png


  • Standard DL1a MC:
- Training:
/fefs/aswg/data/mc/DL1/AllSky/20220511_allsky_std/TrainingDataset/dec_4822/GammaDiffuse/
/fefs/aswg/data/mc/DL1/AllSky/20220511_allsky_std/TrainingDataset/dec_4822/Protons/
- Testing (Ring-like MC, offset 0.4 deg):
/fefs/aswg/data/mc/DL1/AllSky/20220511_allsky_std/TestingDataset
- Testing (Diffuse MC, only four nodes close to the source path in the sky):
/fefs/aswg/data/mc/DL1/AllSky/20220527_src2_diffgamma/TestingDataset/GammaDiffuse/dl1_20220527_src2_diffgamma_{NODE}_merged.h5


  • Tuned DL1b MC:
-NSB tuning parameters:
"image_modifier": {
     "increase_nsb": true,
     "extra_noise_in_dim_pixels": 0.919,
     "extra_bias_in_dim_pixels": 0.298,
     "transition_charge": 8,
     "extra_noise_in_bright_pixels": 0.972,
     "increase_psf": false,
     "smeared_light_fraction": 0.0
   },


- Training:
/fefs/aswg/data/mc/DL1/AllSky/20220524_dec_4822_tuned/TrainingDataset/dec_4822/GammaDiffuse/dl1_20220524_dec_4822_tuned_dec_4822_GammaDiffuse_merged.h5
/fefs/aswg/data/mc/DL1/AllSky/20220524_dec_4822_tuned/TrainingDataset/dec_4822/Protons/dl1_20220524_dec_4822_tuned_dec_4822_Protons_merged.h5
- Testing:
  • gamma point-like
/fefs/aswg/data/mc/DL1/AllSky/20220524_dec_4822_tuned/TestingDataset/dl1_20220524_dec_4822_tuned_node_theta_{NODE}__merged.h5
  • gamma diffuse
For the diffuse-gamma test sample, were used DL1 tuned with slightly different values for the nsb parameters:
"image_modifier": {
  "increase_nsb": true,
  "extra_noise_in_dim_pixels": 0.979,
  "extra_bias_in_dim_pixels": 0.339,
  "transition_charge": 8,
  "extra_noise_in_bright_pixels": 1.133,
  "increase_psf": false,
  "smeared_light_fraction": 0.0

}

/fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/data/mc/ALLsky/TestDiffuse/DL1/dl1_20220527_src2_diffgamma_node_theta_{NODE}__merged.h5
  • Random Forests (source independent):
- Models trained in lstchain v0.9.6, i.e. with "az_tel", "alt_tel" RF features
- cfg file: /fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/lstchain_config_2022-05-24.json
 /fefs/aswg/data/models/AllSky/20220524_dec_4822_tuned/dec_4822
  • DL2 MC (testing only for IRFs):
- gamma point-like
/fefs/aswg/data/mc/DL2/AllSky/20220524_dec_4822_tuned/TestingDataset/dec_4822/{NODE}/dl2_20220524_dec_4822_tuned_node_theta_{NODE}__merged.h5
- gamma diffuse
/fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/data/mc/ALLsky/TestDiffuse/DL2/dl2_20220527_src2_diffgamma_node_theta_{NODE}__merged.h5

DL1 data

  • DL1a files produced by LSTOSA (lstchain v0.9)
  • lstchain v0.9 tailcut8-4 (with cleaning based on pedestal RMS, dynamical cleaning)
  • cfg /fefs/aswg/data/real/DL1/{date}/v0.9/tailcut84/log/lstchain_config_tailcut84_v092.json
/fefs/aswg/data/real/DL1/{date}/v0.9/tailcut84/

DL2 data

  • whole data sample before data selection
/fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/data/DL2_ALLsky

DL3 data

  • Fixed cuts:
"intensity": [50, Infinity],
"width": [0, Infinity],
"length": [0, Infinity],
"r": [0, 1],
"wl": [0.1, 1],
"leakage_intensity_width_2": [0, 1.0],
"event_type": [32, 32]

IRFs

  • Full enclosure
  • Four different IRFs merged with DL2 data depending on the run zenith angle
  • Global cuts optimized on Crab data used
  • Produced in lstchain v0.9.6
  • IRF path: /fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/data/mc/ALLsky/TestDiffuse/IRFs
irf_20220527_src2_diffgamma_node_theta_{node}__merged.fits.gz

IRFs Aeff.png

DL3 files

/fefs/aswg/workspace/giorgio.pirola/LST_analysis/lhaaso_pipe/data/DL3_ALLsky/DiffusTestGammas

High-level analysis

Theta2 distribution

- Energy dependent gammaness cuts optimized on Crab detection significance + global theta2 cut (also optimized on Crab significance, but not bin-wise)
- 3 OFF regions
- Three energy bins (blind), for E>3TeV we have 3.54 sigma

Theta2 fromDL2.png

Excess plot fromDL2.png

1D Spectral analysis

  • 2 hypothesis:
- point-like source: theta<0.2deg (Cut optimized on Crab data)
- extended source: theta<0.26deg (UL reported in the LHAASO paper)
  • Performed with gammapy-v0.19
bkg_maker = ReflectedRegionsBackgroundMaker(exclusion_mask=exclusion_mask,
                                           min_distance=1 * u.rad, # Minimum distance from input region
                                           max_region_number=2 # Maximum number of OFF regions
                                          )
safe_mask_masker = SafeMaskMaker(methods=["aeff-max"], aeff_percent=10)
e_reco_min = 0.1
e_reco_max = 100
e_true_min = 0.01
e_true_max = 100
energy_axis = MapAxis.from_energy_bounds(
   e_reco_min, e_reco_max, 
   nbin=2, per_decade=True, 
   unit="TeV", name="energy"
)
energy_axis_true = MapAxis.from_energy_bounds(
   e_true_min, e_true_max, 
   nbin=5, per_decade=True, 
   unit="TeV", name="energy_true"
)

  • Spectral fitting of stacked LST dataset
  • Performed for point-like source assumption
  • LST-1 data alone: Power-law spectral model
  • Joint likelihood forward folding with LHAASO SED data points: Exponential cutoff power-law spectral model

Power-law model of LST-1 data

Point-like Table powerlaw point allsky.png

ALLsky vs old SEDs point.png


Extended Table powerlaw extended allsky.png

ALLsky vs old SEDs extended.png

Joint likelihood forward folding with LHAASO flux points

Point-like Table ECPL point allsky.png

ALLsky vs old SEDs joint point.png


Extended Ttable ECPL extended allsky.png

ALLsky vs old SEDs joint extended.png

Comparison with Jakub's results

Point-like

ALLsky giorgio vs Jakub SEDs.png

ALLsky giorgio vs Jakub SEDs joint.png

Analysis of systematic

2D Analysis

For the bkg modelling, it has been used a package developed for creati radial acceptance model to be used for IACT analysis with gammapy: https://github.com/mdebony/acceptance_modelisation

  • Acceptance calculation

The acceptance evaluation has been performed in 2 ways:

- Full range
- Dividing the sample in zenith bins:
zenith_bins = np.array([[0, 23.5],
                      [23.5, 30.0],
                      [30.0, 37.6],
                      [37.6, 47.8],
                      [47.8, 55]]
                     )#deg


In both cases, an exclusion mask was used on the 3 sources identified in the ROI:

LHAASO_Ra = 317.22 * u.deg
LHAASO_DEC = 51.95 * u.deg
LHAASO_source = SkyCoord(LHAASO_Ra, LHAASO_DEC, frame='icrs')
RA_4GL=317.01670*u.deg
DEC_4GL=51.9276*u.deg
source_4GL=SkyCoord(RA_4GL, DEC_4GL, frame='icrs')
RA_HARD=317.59578276*u.deg
DEC_HARD=51.79427743*u.deg
source_HARD=SkyCoord(RA_HARD, DEC_HARD, frame='icrs') 
exclude_regions=[CircleSkyRegion(center=LHAASO_source,
                                radius=0.2*u.deg),CircleSkyRegion(center=source_4GL,
                                radius=0.2*u.deg),CircleSkyRegion(center=source_HARD,
                                radius=0.2*u.deg)]

Exclusion region.png

# Define the binning of the model
e_min, e_max = 3*u.TeV, 100*u.TeV
size_fov = 2.5*u.deg
offset_axis_acceptance = MapAxis.from_bounds(0.*u.deg, size_fov, nbin=6, name='offset')
energy_axis_acceptance = MapAxis.from_energy_bounds(e_min, e_max, nbin=10, name='energy')
acceptance_model_creator = RadialAcceptanceMapCreator(energy_axis_acceptance, 
                                                     offset_axis_acceptance,
                                                     exclude_regions=exclude_regions,
                                                     oversample_map=100)
acceptance_model = acceptance_model_creator.create_radial_acceptance_map(obs_collection)

Acceptance.png

  • Skymaps

The skymaps have been produce for the low and the high energy-range, for both the whole-sample and the zenith-wise acceptance cases. It were also noticed few differences if using or not a safe mask:

config.datasets.safe_mask.methods=["aeff-max"] 
config.datasets.safe_mask.parameters={'aeff_percent':10}
  • Acceptance model built using the full dataset:
- Without safe mask:

Skymaps full no mask.png

File:Sig dist full no mask.png

- With safe mask:

Skymaps full yes mask.png

File:Sig dist full yes mask.png

  • Acceptance model built dividing the dataset Zenith-wise:
- Without safe mask:

Skymaps zd no mask.png

Sig dist zd no mask.png

- With safe mask:

Skymaps zd yes mask.png

Sig dist zd yes mask.png