Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Nckw add example teststat plots #845

Merged
merged 5 commits into from
Jul 6, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 13 additions & 2 deletions docs/part3/commonstatsmethods.md
Original file line number Diff line number Diff line change
Expand Up @@ -599,10 +599,21 @@ The distributions of the test-statistic can also be plotted, at each value in th
python test/plotTestStatCLs.py --input mygrid.root --poi r --val all --mass MASS
```

The resulting output file will contain a canvas showing the distribution of the test statistic background only and signal+background hypothesis at each value of **r**. Use `--help` to see more options for this script.
The resulting output file will contain a canvas showing the distribution of the test statistic background only and signal+background hypothesis at each value of **r**. Use `--help` to see more options for this script.

!!! info
If you used the TEV or LEP style test statistic (using the commands as described above), then you should include the option `--doublesided`, which will also take care of defining the correct integrals for $p_{\mu}$ and $p_{b}$.
If you used the TEV or LEP style test statistic (using the commands as described above), then you should include the option `--doublesided`, which will also take care of defining the correct integrals for $p_{\mu}$ and $p_{b}$. Click on the examples below to see what a typical output of this plotting tool will look like when using the LHC test statistic, or TEV test statistic.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not necessarily for this PR, but I wonder if it is worth making this behaviour default, based on the test stat that is being used, rather than requiring the user to do it manually?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The reason I would leave this option open is that we did have a workflow that was specifically for spin Hypotheses and iirc it was not quite the same as using the TeV test stat but was double sided. Honestly, if you don't use the option when you should, you'll see pretty quickly from the resulting figure that you should have used it. In any case, this script doesn't know about the test stat but I think what you're saying is you we could pass the test-stat in as an option and have the plot figure it out?


<details>
<summary><b>qLHC test stat example</b></summary>
![](images/exampleLHC.jpg)
</details>

<details>
<summary><b>qTEV test stat example</b></summary>
![](images/exampleTEV.jpg)
</details>


## Computing Significances with toys

Expand Down
Binary file added docs/part3/images/exampleLHC.jpg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added docs/part3/images/exampleTEV.jpg
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
20 changes: 20 additions & 0 deletions test/plotTestStatCLs.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
ROOT.gROOT.SetBatch(1)
ROOT.gROOT.ProcessLine(".L $CMSSW_BASE/src/HiggsAnalysis/CombinedLimit/test/plotting/hypoTestResultTree.cxx")
ROOT.gROOT.ProcessLine(".L $CMSSW_BASE/src/HiggsAnalysis/CombinedLimit/test/plotting/qmuPlot.cxx")
ROOT.gStyle.SetOptStat(0)

from ROOT import hypoTestResultTree, q0Plot, qmuPlot

Expand Down Expand Up @@ -42,6 +43,14 @@
type="str",
help="poi values, comma separated (type all to make a plot for every value found in the file)",
)
parser.add_option(
"",
"--sub-label",
default="",
type="str",
dest="sub_label",
help="change sub-label from q to q_sub_label (doesn't apply if using --signif option",
)
parser.add_option(
"-r",
"--rebin",
Expand Down Expand Up @@ -99,6 +108,14 @@
default=False,
help="If --signif, Force plotting of both distributions (including for signal injected).",
)
parser.add_option(
"",
"--save-as-pdf",
dest="save_as_pdf",
action="store_true",
default=False,
help="Save plots as pdfs as well as Canvases in root files.",
)
(options, args) = parser.parse_args()

if options.quantileExpected >= 0:
Expand Down Expand Up @@ -167,8 +184,11 @@ def findPOIvals(fi, m):
options.rebin,
options.expected,
options.quantileExpected,
options.sub_label,
)
canvases.append(can.Clone())
if options.save_as_pdf:
can.SaveAs(options.output + "_" + can.GetName() + ".pdf")
ft.Close()
ofile = ROOT.TFile(options.output, "RECREATE")
for can in canvases:
Expand Down
19 changes: 16 additions & 3 deletions test/plotting/qmuPlot.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ TCanvas *q0Plot(float mass, std::string poinam , float poival, int rebin=0, bool



TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int invert=0,int rebin=0, int runExpected_=0, double quantileExpected_=0.5) {
TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int invert=0,int rebin=0, int runExpected_=0, double quantileExpected_=0.5,std::string testStatSublabel="") {
if (gFile == 0) { std::cerr << "You must have a file open " << std::endl; return 0; }
TTree *t = (TTree *) gFile->Get("q");
if (t == 0) { std::cerr << "File " << gFile->GetName() << " does not contain a tree called 'q'" << std::endl; return 0; }
Expand All @@ -213,6 +213,16 @@ TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int
TH1F *qB;
TH1F *qS;

// check if should use --doublesided (mode==1) option
if (mode==0) {
TH1F *qTest;
t->Draw("2*q>>qTest","weight*(type==-1)");
qTest = (TH1F*) gROOT->FindObject("qTest")->Clone();
if (qTest->Integral(1,qTest->FindBin(0)) > 0.3*qTest->Integral()) {
std::cout << "WARNING -- It looks like you are using either the TEV or LEP style test-statistic. If so, you should use the option --doublesided" << std::endl;
}
}

if (mode==0) t->Draw("max(2*q,0)>>qB","weight*(type==-1)");
else t->Draw("2*q>>qB","weight*(type==-1)");

Expand Down Expand Up @@ -321,7 +331,9 @@ TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int
}
leg1->AddEntry(qO, "observed value", "L");

TLegend *leg2 = new TLegend(.63,.67,.93,.48);
TLegend *leg2;
if (mode==0) leg2 = new TLegend(.63,.67,.93,.48);
else leg2 = new TLegend(.13,.87,.43,.68);
leg2->SetFillColor(0);
leg2->SetShadowColor(0);
leg2->SetTextFont(42);
Expand Down Expand Up @@ -354,7 +366,8 @@ TCanvas *qmuPlot(float mass, std::string poinam, double poival, int mode=0, int
leg2->Draw();
qB->SetTitle("");
qB->GetYaxis()->SetTitle("");
qB->GetXaxis()->SetTitle(Form("q_{%s}(%s = %g, m_{H} = %g GeV)",poinam.c_str(),poinam.c_str(),poival,mass));
if (testStatSublabel.length()>0) qB->GetXaxis()->SetTitle(Form("q_{%s}(%s = %g, m_{H} = %g GeV)",testStatSublabel.c_str(),poinam.c_str(),poival,mass));
else qB->GetXaxis()->SetTitle(Form("q_{%s}(%s = %g, m_{H} = %g GeV)",poinam.c_str(),poinam.c_str(),poival,mass));
qB->GetXaxis()->SetTitleOffset(1.05);


Expand Down