CLs upper limits using RooStats

Hi,

I hope that this is the right place for my question.

I am trying to set 95% upper limits using CLs on the signal events for the event counting experiment.
I would like to know if the following steps are correct:

  1. Create the example.root file:

[code]void makeChannel0()
{
TFile* example = new TFile(“data/example.root”,“RECREATE”);

TH1F* data = new TH1F(“data”,“data”, 1,0,1);

TH1F* signal = new TH1F(“signal”,“signal histogram (pb)”, 1,0,1);

TH1F* background = new TH1F(“background”,“background histogram (pb)”, 1,0,1);

data->Fill(0.5,69.00);
background->Fill(0.5,74.61);
signal->Fill(0.5,66.00);

example->Write();
example->Close();
}
[/code]

  1. config/example.xml file:

[code]

./config/example_channel0.xml

norm Lumi [/code]
  1. config/example_channel0.xml file:

[code]

<Sample Name="signal" HistoPath=""  HistoName="signal">
  <NormFactor Name="norm" Val="1" Low="0." High="2."  />
</Sample>

<Sample Name="background" HistoPath="" NormalizeByTheory="True" HistoName="background">
  <OverallSys Name="syst1" Low="0.93" High="1.05"/>
</Sample>
[/code]
  1. Create the ModelConfig:
  1. Compute limits:

Could you please tell me if the above steps are correct?

Thanks a lot,
Archil

Hi,

What you are doing seems correct to me. Are you not getting results which do make sense ? If this is the case, please attach the workspace and/or your input xml and root file, so I can investigate the problem, if there is one

Best Regards
Lorenzo

Hi Lorenzo,
Thank you very much for the answer.

I am glad if it seems correct what I am doing :slight_smile:

I do not have a big experience in RooStats and I just wanted to be sure that it is correct what I am doing…

Also I would like to ask you one question:
If I include too many systematic uncertainties sources in the ModelConfig, than the results using Frequentist and Asymptotic calculator are significantly different. As I understund it is expected. Could you confirm that this difference is expected?

Best Regards,
Archil

Hi,

The result should not be too different, especially in case you have a model where the statistics is low (and the asymptotic limits are valid). If you observe in that case a big difference i maybe a problem in the way the model is defined. For example, have you correctly defined in the ModelConfig the nuisance parameters ? Do you have also global observables ? They need to be defined in order to be randomised in the FrequentistCalculator

Best Regards

Lorenzo

Hi Lorenzo,
Thanks a lot for helping me.

Yes, I have the global observables and nuisance parameters in ModelConfig. This is my ModelConfig with few nuisance parameters:

=== Using the following for ModelConfig === Observables: RooArgSet:: = (obs_x_channel0) Parameters of Interest: RooArgSet:: = (norm) Nuisance Parameters: RooArgSet:: = (Lumi,alpha_BJesUnc,alpha_btag,alpha_jeff,alpha_jer,alpha_jvf,alpha_mistag) Global Observables: RooArgSet:: = (nominalLumi,nom_alpha_BJesUnc,nom_alpha_btag,nom_alpha_jeff,nom_alpha_jer,nom_alpha_jvf,nom_alpha_mistag) PDF: RooProdPdf::model_channel0[ lumiConstraint * alpha_BJesUncConstraint * alpha_btagConstraint * alpha_jeffConstraint * alpha_jerConstraint * alpha_jvfConstraint * alpha_mistagConstraint * alpha_A_NormConstraint * alpha_B_NormConstraint * channel0_model(obs_x_channel0) ] = 74.0024

If I include few systematic uncertainties sources in the ModelConfig, then expected and observed limits using Frequentist and Asymptotic calculator are similar.
But if I include too many systematic uncertainties sources, then observed limit using Frequentist calculator is much higer than the one using Asymptotic calculator, but expected limits are similar using both clalculator.

Here are the files which I use in case of few systematics:

  1. file for creating example.root file:

[code]void makeChannel0()
{
TFile* example = new TFile(“data/example.root”,“RECREATE”);

TH1F* data = new TH1F(“data”,“data”, 1,0,1);
TH1F* signal = new TH1F(“signal”,“signal histogram (pb)”, 1,0,1);
TH1F* background_A = new TH1F(“background_A”,“background histogram (pb)”, 1,0,1);
TH1F* background_B = new TH1F(“background_B”,“background histogram (pb)”, 1,0,1);
TH1F* background_C = new TH1F(“background_C”,“background histogram (pb)”, 1,0,1);
TH1F* background_D = new TH1F(“background_D”,“background histogram (pb)”, 1,0,1);
TH1F* background_E = new TH1F(“background_E”,“background histogram (pb)”, 1,0,1);
TH1F* background_F = new TH1F(“background_F”,“background histogram (pb)”, 1,0,1);
TH1F* background_G = new TH1F(“background_G”,“background histogram (pb)”, 1,0,1);
TH1F* background_H = new TH1F(“background_H”,“background histogram (pb)”, 1,0,1);
TH1F* background_I = new TH1F(“background_I”,“background histogram (pb)”, 1,0,1);

data->Fill(0.5,69.0000);
background_A->Fill(0.5,16.8921);
background_B->Fill(0.5,40.4409);
background_C->Fill(0.5,6.5852);
background_D->Fill(0.5,1.9575);
background_E->Fill(0.5,3.8127);
background_F->Fill(0.5,4.8727);
background_G->Fill(0.5,0.1594);
background_H->Fill(0.5,0.0422);
background_I->Fill(0.5,1.0451);
signal->Fill(0.5,166.91);

example->Write();
example->Close();
}
[/code]

  1. config/example.xml - file:

[code]

./config/example_channel0.xml

norm Lumi alpha_A_Norm alpha_B_Norm [/code] 3. config/example_channel0.xml file: [code]
<Sample Name="background_A" HistoPath="" NormalizeByTheory="True" HistoName="background_A">
  <OverallSys Name="BJesUnc" Low="1.0061" High="1.0642"/>
  <OverallSys Name="btag" Low="0.9694" High="1.0657"/>
  <OverallSys Name="jeff" Low="0.9753" High="1.0247"/>
  <OverallSys Name="jer" Low="0.9413" High="1.0587"/>
  <OverallSys Name="jvf" Low="1.0226" High="1.0228"/>
  <OverallSys Name="mistag" Low="0.9823" High="1.0924"/>
  <OverallSys Name="A_Norm" Low="0.8590" High="1.1410"/>
</Sample>

<Sample Name="background_B" HistoPath="" NormalizeByTheory="True" HistoName="background_B">
  <OverallSys Name="BJesUnc" Low="1.0450" High="1.0255"/>
  <OverallSys Name="btag" Low="1.0690" High="1.0516"/>
  <OverallSys Name="jeff" Low="0.9314" High="1.0686"/>
  <OverallSys Name="jer" Low="0.8250" High="1.1750"/>
  <OverallSys Name="jvf" Low="1.0393" High="0.9988"/>
  <OverallSys Name="mistag" Low="1.0726" High="1.0897"/>
  <OverallSys Name="B_Norm" Low="0.8650" High="1.1350"/>
</Sample>

<Sample Name="background_C" HistoPath="" NormalizeByTheory="True" HistoName="background_C">
  <OverallSys Name="BJesUnc" Low="1.0731" High="1.0983"/>
  <OverallSys Name="btag" Low="1.1072" High="0.8890"/>
  <OverallSys Name="jeff" Low="0.9568" High="1.0432"/>
  <OverallSys Name="jer" Low="0.8769" High="1.1231"/>
  <OverallSys Name="jvf" Low="1.0463" High="1.0277"/>
  <OverallSys Name="mistag" Low="1.0002" High="1.0524"/>
</Sample>

<Sample Name="background_D" HistoPath="" NormalizeByTheory="True" HistoName="background_D">
 <OverallSys Name="BJesUnc" Low="1.0031" High="1.0001"/>
 <OverallSys Name="btag" Low="1.1350" High="0.8650"/>
 <OverallSys Name="jeff" Low="1.0000" High="1.0000"/>
 <OverallSys Name="jer" Low="0.3510" High="1.6490"/>
 <OverallSys Name="jvf" Low="1.0000" High="1.0000"/>
 <OverallSys Name="mistag" Low="0.9739" High="1.0258"/>
</Sample>

<Sample Name="background_E" HistoPath="" NormalizeByTheory="True" HistoName="background_E">
 <OverallSys Name="BJesUnc" Low="0.9755" High="0.9928"/>
 <OverallSys Name="btag" Low="0.9671" High="0.9279"/>
 <OverallSys Name="jeff" Low="0.9482" High="1.0518"/>
 <OverallSys Name="jer" Low="0.8532" High="1.1468"/>
 <OverallSys Name="jvf" Low="0.9488" High="0.9982"/>
 <OverallSys Name="mistag" Low="0.9213" High="1.0132"/>
</Sample>

<Sample Name="background_F" HistoPath="" NormalizeByTheory="True" HistoName="background_F">
 <OverallSys Name="BJesUnc" Low="0.9792" High="0.9869"/>
 <OverallSys Name="btag" Low="1.0322" High="0.9288"/>
 <OverallSys Name="jeff" Low="0.9799" High="1.0201"/>
 <OverallSys Name="jer" Low="0.9873" High="1.0127"/>
 <OverallSys Name="jvf" Low="0.9938" High="0.9794"/>
 <OverallSys Name="mistag" Low="0.9895" High="0.9707"/>
</Sample>

<Sample Name="background_G" HistoPath="" NormalizeByTheory="True" HistoName="background_G">
 <OverallSys Name="BJesUnc" Low="1.0019" High="1.0006"/>
 <OverallSys Name="btag" Low="1.0620" High="0.9243"/>
 <OverallSys Name="jeff" Low="1.0000" High="1.0000"/>
 <OverallSys Name="jer" Low="0.6301" High="1.3699"/>
 <OverallSys Name="jvf" Low="1.0000" High="1.0000"/>
 <OverallSys Name="mistag" Low="1.0069" High="0.9931"/>
</Sample>

<Sample Name="background_H" HistoPath="" NormalizeByTheory="True" HistoName="background_H">
 <OverallSys Name="BJesUnc" Low="1.0025" High="1.0161"/>
 <OverallSys Name="btag" Low="1.1241" High="0.8483"/>
 <OverallSys Name="jeff" Low="0.9566" High="1.0434"/>
 <OverallSys Name="jer" Low="0.9130" High="1.0870"/>
 <OverallSys Name="jvf" Low="1.0678" High="1.0073"/>
 <OverallSys Name="mistag" Low="0.9905" High="0.9859"/>
</Sample>

<Sample Name="background_I" HistoPath="" NormalizeByTheory="True" HistoName="background_I">
 <OverallSys Name="BJesUnc" Low="1.0103" High="1.0538"/>
 <OverallSys Name="btag" Low="0.9790" High="1.0465"/>
 <OverallSys Name="jeff" Low="0.9488" High="1.0512"/>
 <OverallSys Name="jer" Low="0.9746" High="1.0254"/>
 <OverallSys Name="jvf" Low="1.0399" High="1.0399"/>
 <OverallSys Name="mistag" Low="1.0298" High="1.0188"/>
</Sample>

[/code]

If I use the files shown above, then I get the follownig results:
Asymptotic calculator:

The computed upper limit is: 0.154325 +/- 0 Expected upper limits, using the B (alternate) model : expected limit (median) 0.184132 expected limit (-1 sig) 0.135239 expected limit (+1 sig) 0.250728 expected limit (-2 sig) 0.0978015 expected limit (+2 sig) 0.300834
Frequentist calculator (10000 toy):

The computed upper limit is: 0.177149 +/- 0.00358786 Expected upper limits, using the B (alternate) model : expected limit (median) 0.183857 expected limit (-1 sig) 0.136362 expected limit (+1 sig) 0.250136 expected limit (-2 sig) 0.112155 expected limit (+2 sig) 0.295055
In this case the results are similar and everythig is fine, but if I use the xml files which are in attachment (all systematic uncertainties sources are included) then I get the following results:
Asymptotic calculator:

The computed upper limit is: 0.18238 +/- 0 Expected upper limits, using the B (alternate) model : expected limit (median) 0.198307 expected limit (-1 sig) 0.145693 expected limit (+1 sig) 0.282874 expected limit (-2 sig) 0.112837 expected limit (+2 sig) 0.34469
Frequentist calculator (10000 toy):

The computed upper limit is: 0.314777 +/- 0.00241345 Expected upper limits, using the B (alternate) model : expected limit (median) 0.207605 expected limit (-1 sig) 0.161386 expected limit (+1 sig) 0.299007 expected limit (-2 sig) 0.0788636 expected limit (+2 sig) 0.39817
Here we can see that observed limits are significantly different.

Sorry that my post is too big :frowning:
but it’s very important for me to understand what’s happening and if the results are correct.

Best regards,
Archil
example.tar.gz (4.29 KB)

Hi,

Can you please post all the files needed to run your model. You can post your workspaces, or a tar file with all the needed files (including the input histogram files), that I can run hist2workspace.
A simple xml file is not enough

Best Regards

Lorenzo

Hi Lorenzo,

I have attached all the files which I use.

when I introduce the systematics, for example this one:

I put the value corresponding to “jvf_down” as Low attribute of OverallSys and the value corresponding to “jvf_up” - as High attribute of OverallSys.
Hope it is a correct way.

Thanks a lot.

Best Regards,
Archil
example.tar.gz (201 KB)

Hi,

Thank you for the files. Yes, I can reproduce now the difference between the asymptotic and the frequentist calculator. It looks that the test statistic distribution obtained from the toys in the case of the signal plus background model are not following a chi-square distribution. It looks like that, when generating toys with the signal, the model is not sensitive to the presence of the signal.
Do you have also the xml data files with less systematic uncertainty which give a result compatible between asymptotic and frequentist calculators ?

Best Regards

Lorenzo

Hi Lorenzo,

Thanks a lot for the reply.

I have attached all the files with less systematic uncertainty.

I think that the systematic uncertainty on background is too big, it’s about 25%, and this could be the reason that the model is not sensitive to the presence of the signal.

Best Regards,
Archil
example_fewSyst.tar.gz (146 KB)

Hi,

Thank you for your new files with few systematics. In this case I see the test statistic distributions for the S+B toys following a chi2 distribution and in did I get a result as expected. I will investigate why we have this effect in case of many systematics.

Lorenzo

Dear Lorenzo,
Thanks a lot for the reply.

Then I will wait the results which you get when investigate why we have this effect.

Sorry for disturb you,

Best Regards,
Archil

Hi,

I have investigate your workspace and found that , when generating toys with a poi value (e.g. “norm” = 0.5), you get a best fit value of norm significantly lower.
It looks like all your toys are biased. I don;t know if this is a problem in the model or the generation. This explains what you observe.
Can you maybe try to remove one by one the systematics to see if this is due to one of them ? Otherwise it is difficult to find out what could bias your toys

Best Regards

Lorenzo

Hi Lorenzo,

I have tried to decrease the number of included systematic sources on the background and I see that if it is included any 6 source of systematics from all of them, then the results using both calculator are close to each other.
If I increase the number of included systematic sources, then the difference also increases between the results using different calculator.

I have attached the files included 6 and 9 systematics sources.

If I use 9 source of systematics but I set every systematics on background to be 5%, then the results using both calculator are close to each other. these files are also attached.

I do not understand the results :frowning:

Thanks a lot.

Best Regards,
Archil
example_5percent.tar.gz (179 KB)
example_9syst.tar.gz (184 KB)
example_6syst.tar.gz (183 KB)

Hi,

Yes I can reproduce the results. I see well the case by looking at the fitted poi parameter of every toys (see the top-left histograms in the attached plots in case of example_5percent and example_9syst).
The toys are generated with POI=0.5 and in the first case the fitted distribution of the POI is in average 0.5 while in the second is clearly distorted and is not a gaussian with mean 0.5 as should be.
Why this is happening I do not know, it is something that needs to be understood

Best Regards

Lorenzo
toys_5percent.pdf (18.2 KB)
toys_9syst.pdf (18.7 KB)

Hi Lorenzo,

Thanks for the plots.
Could you tell me how i can obtain the same plots?

Yes now I see the difference.
probably problem is in the generation of the global observables, but if the model is same, why we get so different disstribution of fitted POI…

Best Regards,
Archil

Hi,

You can make these plots using the detailed output capabilities of the RooStats test statistics. The detail output will contained the fitted parameter values for each toy and the best fit likelihood values. I have committed a new version of StandardHypoTestInvDemo (root.cern.ch/gitweb?p=root.git;a … 861d0874df ) with this capabilities.
You have to do:

.L StandardHypoTestInvDemo.C
enableDetailedOutput=true;
StandardHypoTestInvDemo()

The macro will create an output root file with the HypoTestInvResult object, which has the detailed output information. I attach a macro that I have used to read that information.

Best Regards

Lorenzo
readHypoTestInvDetailedOutput.C (3.09 KB)

Hi Lorenzo,

Thanks a lot for the instructions and the macro.

I think that in case of many asymmetric systematics sources the mean value of the distribution of should not be equal to POI, because in the fit of likelihood the nuisence parameters are constrained: .
from the following plots we can see that if I include asymetric systematic sources then the mean value of the distribution of is not equal to POI=0.5:

the files which i used are in attachment.

Concerning the results using Asimptotic calculator: I am not sure but I think that this calculator uses the following asymptotic formulae [1] :

but this is the special case of the formulae:

We get (56) from (55) if . If we do not have this special case then using of the formula (56) is not correct.

is the mean value of the distribution of . In case of few systematic sources we have the special case and the results using both calculator are similar.

Could you please tell me which formula is used by Asymptotic Calculator?

Best Regards,
Archil
[1] arxiv.org/pdf/1007.1727v3.pdf
asymmetric_syst.tar.gz (193 KB)

Hi ,

The asymptotic calculator uses in your case the formula defined in paragraph 3.7 of the above cited paper.
All the asymptotic formulae are based on the Wald approximation (i.e. asymptotically all the distributions of the parameters are gaussians). It is clear the approximation is not valid if your model is biased as in your case.

I think a small asymmetry in the systematic error is not a problem, but in your case, looking at the xml files I see
something which does not make sense:

<Sample Name="background_A" HistoPath="" NormalizeByTheory="True" HistoName="background_A">

The Low values are larger than 1. This I think makes your systematic effect to bias the model. Is this really wanted ?

For example, whatever you change that particular systematic sources (e.g el_trigSF ), you will have always more background. It is clear that this model is biased !
One should avoid bias in the models, otherwise you would need to correct for it at the end.
It is also clear the asymptotic calculator does not work in a biased model.

Best Regards

Lorenzo

Hi Lorenzo,

Thanks a lot, now I understand what is the problem.
I have recalculated the systematics on total background:

<OverallSys Name="btag" Low="1.03722" High="1.01615"/> <OverallSys Name="eer" Low="1.02959" High="1.01869"/> <OverallSys Name="ees" Low="1.01571" High="1.07095"/> <OverallSys Name="el_trigSF" Low="1.08758" High="1.06476"/> <OverallSys Name="el_idSF" Low="1.0249" High="1.07381"/> <OverallSys Name="el_recSF" Low="1.03171" High="1.08187"/> <OverallSys Name="letp_trigSF" Low="1.01431" High="1.04232"/> <OverallSys Name="flavor_comp" Low="0.982308" High="1.10519"/> <OverallSys Name="flavor_response" Low="0.990957" High="1.09899"/>
and now I can clearly see that I have more background in any case (most of the systematic sources always have a variation in positive direction).

Thank you very much again, conversation with you was very very productive and helpful for me.

and sorry for disturb you.

Best Regards,
Archil