Energy-dependent pull correction

[sathanas31]

When applying an energy-dependent pull correction (eg median_1d), you weight the MC to get the (weighted) median pull. If you are using the power_law energy pdf, you need to provide the gamma in the e_pdf_dict which is then used in the weight_mc().

As it is right now, the pull_dict is created by taking the llh_dict["llh_energy_pdf"], in which you need to pass the gamma in addition to the energy_pdf_name when setting up the mh_dict so that the pulls can be calculated. Alternatively, you can use the inj_dict["injection_energy_pdf"] when creating the pull_dict in the add_angular_error_modifier()

BaseAngularErrorModifier.create(
            season,
            self.inj_dict["injection_energy_pdf"],
            self.floor_name,
            self.pull_name,
            gamma_precision=self.llh_dict.get("gamma_precision", "flarestack"),
        )

Neither of the two solutions is what we need, so the question is do we need to weight the MC to get the pulls? If yes, then the pull_dict needs to be created in such a way that the gamma is provided, but without requiring that the energy pdf is the same for both llh and injection. Any ideas on this?

To Reproduce
Run any analysis with llh_dict["pull_name"] = "median_1d"
Either pass the gamma in the llh_dict:

llh_energy = {"energy_pdf_name": "power_law",
                  "gamma": gamma}
llh_dict["llh_energy_pdf"] = llh_energy

or change the BaseAngularErrorModifier as above

Expected behavior
It will run, after some errors are fixed :P

Additional context
Errors to fix so it can run

1. in get_mc method: floor_dict["season"].mc_path
2. when reading the pull pickle in the BaseMedianAngularErrorModifier: open(self.pull_name, "rb")

[mlincett]

I am not sure I grasp all the details. I somewhat understand the expected behaviour but I don't quite understand the current behaviour and why it is problematic.

Maybe you can elaborate about the "neither of the two solutions is what we need".

In general, in am in favour in tidying up the configuration dictionaries in any sensible way :)

[sathanas31]

The problem lies in how the BaseAngularErrorModifier is called in the minimizer, namely the e_pdf_dict that goes into the pull_dict is taken from the llh_dict["llh_energy_pdf"] that you pass when you create the MinimisationHandler.

Now if you choose the power_law as the energy pdf, then you need the gamma in the e_pdf_dict in order to weight the MC, which is used to calculate the median pull (eg see here). However, you don't pass the gamma in the llh_dict["llh_energy_pdf"] since you want to fit it.

The alternative method, where the inj_dict["injection_energy_pdf"] is taken instead for the e_pdf_dict of the BaseAngularErrorModifier, would work because the injection_energy_pdf has the gamma in the dictionary. However, this restricts you to use the power_law also for injection, so you won't be able to inject with a different spectrum. This configuration works for me because I inject with a power law as well, but it's not preferrable in a general sense due to this restriction.

Another thing we could do, is create the pull_dict separately without taking the llh_energy_dict as it is in order to accomodate for when power_law is used as the llh energy pdf. In any case, my main question is more on the conceptual side of things: do we need to weight the MC to get the (energy-dependent) pulls? If weighting isn't needed, then the problem is solved by itself since weight_mc() is the one that requires the gamma

[mlincett]

Thank you for the explanation, I better understand the context and use case now.

My understanding of the pull correction is that it is evaluated by comparing the true directional reco error with the angular error estimator. In this sense, you need to weight the MC because, similarly to other cases reco vs true , the resulting PDFs will depend on the assumed spectrum (since the PSF is energy-dependent).

I believe the reason the 1D pull correction requires a fixed gamma in the likelihood is that, indeed, the pull (and hence the likelihood) depends on the gamma. To consistently fit gamma, one would need to re-evaluate the pull at every point in the gamma space, which I believe is what the 2D pull feature is there for.

So my question is: do you actually want to fit a free gamma parameter when your likelihood embeds a component evaluated with a fixed gamma? Are there specific reasons for not using a 2D correction?

If the answer to these questions is actually yes, then I suppose BaseAngularErrorModifier.create does not really care about whether you pass a inj_dict or llh_dict to it, as long as it contains the necessary fields.

[sathanas31]

The reason I chose the fixed-gamma pull correction is because I'm currently testing only the scenario of gamma = 2 for a subsample of catalog sources, so that I can have a direct comparison with the sensitivities I get with the KDEs (KDEs inherently account for what the pull correction does). I agree that in any other case where you're fitting the gamma the 2D energy-dependent correction is appropriate (eg the median_1d_e BaseAngularErrorModifier).

Nevertheless, if everyone else agrees that weighting the MC is the sensible thing to do, I would suggest that the e_pdf_dict that goes into the pull_dict be set separately in the minimizer

add_angular_error_modifier(self, season):
        e_pdf_dict  = {"energy_pdf_name": self.llh_dict["llh_energy_pdf"],
                               "gamma": self.inj_dict["injection_energy_pdf"]["gamma"],
                             }
        return BaseAngularErrorModifier.create(
            season,
            e_pdf_dict,
            self.floor_name,
            self.pull_name,
            gamma_precision=self.llh_dict.get("gamma_precision", "flarestack"),
        )

in order to be able to inject with different spectrum. This change would be applicable to the fixed-gamma, energy-dependent pull correction, the energy- and declination-dependent pull correction, as well as all the quantile_floor_0d

@JannisNe any comments?