Extracting systematics from profile

Roofit and RooStats
1

dear experts,
I discussed this with several people, but I would need a confirmation that what i am doing is consistent enough.
Basically, I am trying to extract some single source systematic uncertainty associated to one measurement. In order to extract the different systematics, I am proceeding as follow:

  1. start with the histogram templates for signal (nominal and systematics) and data driven background and produce the Roofit workspace using the Histfactory machinery.

  2. fix all nuisance parameters to zero, fit and extract the 68% confidence interval.
    -> This is the nominal signal strength factor, that i label: mu_nom_c, mu_nom_u, mu_nom_d (central, upper, lower)

  3. leave all nuisance parameters free within [-2,2], perform the fit and extract again the associated 68% confidence interval.
    -> I label these: mu_sys_c, mu_sys_u, mu_sys_d (central, upper, lower)
    -> I have the set of associated values for the nuisance parameters.

  4. loop over all nuisance parameters, for each of them, fix its value to the value found in step 3) and perform the fit w.r.t. all other ones left free and extract the associated 68% CL interval. I call this:
    -> I label these: mu_alphaI_c, mu_alphaI_u, mu_alphaI_d (central, upper, lower)

  5. fix all nuisance parameters to their value found in step 3) and perform the fit -> this allows me to extract the stat only uncertainty: mu_stat_c, mu_stat_u, mu_stat_d

then i compute the systematics (total and single sources ) as follow:
a) total systematics:

tot_sys_u = ( |(mu_sys_c - mu_sys_u)^2 - (mu_stat_c - mu_stat_u)^2 |)^(1/2)
tot_sys_d = ( |(mu_sys_c - mu_sys_d)^2 - (mu_stat_c - mu_stat_d)^2 |)^(1/2)

b) Luminosity:

lumi_sys_u = ( |(mu_lumi_c - mu_lumi_u)^2 - (mu_stat_c - mu_stat_u)^2 |)^(1/2)
lumi_sys_d = ( |(mu_lumi_c - mu_lumi_d)^2 - (mu_stat_c - mu_stat_d)^2 |)^(1/2)

c) Nuisance parameter associated systematics:

alpha_sys_u = ( | (mu_sys_c - mu_sys_u)^2 - (mu_alpha_c  - mu_alpha_u)^2 |)^(1/2);
alpha_sys_d = ( | (mu_sys_c - mu_sys_d)^2 - (mu_alpha_c  - mu_alpha_d)^2 |)^(1/2);

To do this on the technical aspect, I am using Minuit2 like:
mc->GetPdf()->fitTo(*data,Save(kTRUE),InitialHesse(true), Hesse(true),Minos(kTRUE), RooFit::Minimizer(“Minuit2”,“minos”), Strategy(1), PrintLevel(3), Save(true) );

and to extract the confidence interval, i use:
ProfileLikelihoodCalculator* PL = new ProfileLikelihoodCalculator(*data,mc);
PL->SetConfidenceLevel(0.68);
LikelihoodInterval interval = PL->GetInterval();
double intervalValue = POI->getVal();
double intervalError = POI->getError();
double intervalLow = interval->LowerLimit(*POI);
double intervalUp = interval->UpperLimit(*POI);

thanks in advance.

2

Hi,

Probably your procedure is correct. It is anyway just an assumption to call systematic error the contributions of all the nuisance.
One small comment, I would do step 4 by fixing all the nuisance parameters except one (i.e. you profile the parameter of interest and one nuisance parameter) and then, when computing the single systematic uncertainty you subtract the statistical contribution obtained in step 5. This should be also faster.

Best Regards

Lorenzo

3

many thanks for the reply Lorenzo. this clarifies things.

4

hi again,
one point which is still not clear to me is the following.
Assuming i want to give the systematics that is extracted for each single nuisance parameter in percent.
I am not sure anymore whether i should give this as:

  1. normalised to 1/mu_stat_c * 100
    or
    2)normalised to 1/mu_nom_c * 100
    I naively think that this should be given as 1).
    As a reminder:
    mu_stat_c is derived when fixing all nuisance parameters to their value found in in step 3)
    mu_nom_c is derived with fixing all nuisance parameters to 0.
    thanks for the clarification.