Line: 1 to 1 | ||||||||||||||||||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| ||||||||||||||||||||||||||||||||||||||||||||||
Line: 6 to 6 | ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Higgs Combination and Properties: Long Exercise - CMSDAS 2014 | |||||||||||||||||||||||||||||||||||||||||||||
> > | Higgs combine tutorial based on CMSDAS2014HiggsCombExercise | |||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Line: 27 to 27 | ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Line: 42 to 42 | ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
Intro I: The Higgs at CMSMeasurements of the Higgs properties in CMS exploit the different production and decay modes expected for the SMH boson. | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | The SMH boson is produced through gluon-fusion (ggH, red), vector-boson fusion (VBF, blue), associated vector-boson production (VH, purple - where the Higgs is radiated off an energetic V=W,Z boson), and in associated production with a top-quark pair (ttH, green): ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | The SMH boson is produced through gluon-fusion (ggH, red), vector-boson fusion (VBF, blue), associated vector-boson production (VH, purple - where the Higgs is radiated off an energetic V=W,Z boson), and in associated production with a top-quark pair (ttH, green): ![]() | |||||||||||||||||||||||||||||||||||||||||||||
As can be seen, the latter three modes have very clear signatures that can be targeted experimentally by, e.g., requiring vector boson candidates in the event or top quark candidates. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 57 to 57 | ||||||||||||||||||||||||||||||||||||||||||||||
At the LHC, given the very large gluon luminosities, the SMH boson is mainly produced through gluon-fusion and, at much lower rates, the more distinctive processes: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
At a mass of mH ~ 125 GeV, the SMH boson decays into almost every possibly mode except into a pair of top quarks: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Here, the decay into a pair of photons is also loop-mediated as in the case of gluon-fusion production, except that in this case, the photon is sensitive not only to the colour-charge but also to the electric-charge: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
While the top and W contributions to the diphoton decay, contributions have also been calculated for other quarks and even the tau lepton. In this way, the diphoton decay branching fraction is sensitive no only to the SMH coupling to those SM particles but as well to BSM particles. In terms of decay rates, although they seem very small for some decays, in practice and in order to observe the ZZ or WW decays one needs to take into account the branching fractions of the Z and W. After that, the picture is a little different: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Finally, while it is somewhat straightforward to identify the decay signatures and experimentally select pure samples of H→γγ or H→ττ, it is much more complicated to distinguish between production modes. In an nutshell, only leptonic decays of Vs in VH production are extremely pure with the present amount of data. Both VBF and tt¨¬„H selections inevitably include some gluon-fusion events where there was one or more jets radiated from the gluon or quark lines. With more data, stricter cuts can be applied and purer production tags developed. | |||||||||||||||||||||||||||||||||||||||||||||
> > | Finally, while it is somewhat straightforward to identify the decay signatures and experimentally select pure samples of H→γγ or H→ττ, it is much more complicated to distinguish between production modes. In an nutshell, only leptonic decays of Vs in VH production are extremely pure with the present amount of data. Both VBF and ttì„H selections inevitably include some gluon-fusion events where there was one or more jets radiated from the gluon or quark lines. With more data, stricter cuts can be applied and purer production tags developed. | |||||||||||||||||||||||||||||||||||||||||||||
But for the reasons above we usually refer to "dijet tag" and not "VBF production" when discussing the experimental selections, since such selections will contains some events other than VBF. In CMS, a large amount of production and decay modes are studied exploring different experimental selections targeting the different production modes. The following table summarises the existing and planned analyses: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
There are several notes to this table:
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
From this rich picture of measurements, one can explore the most important aspects of the scalar coupling structure of the new boson. | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | But not all ends with the measurements of 堦BR. | |||||||||||||||||||||||||||||||||||||||||||||
> > | But not all ends with the measurements of σ×BR. | |||||||||||||||||||||||||||||||||||||||||||||
In fact there is a much richer and complex structure in terms of the tensor couplings. This is a more difficult arena where the spin (J) and charge-parity (CP) properties of the boson are involved. | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | The only channel where all information about the boson is available in the detector is the H→ZZ→4e¡° where the leptons are electrons or muons (and not neutrinos). This does not mean that other channels cannot contribute to this picture (like H→WW→e¡°νe¡°ν or H→γγ), but the fact that all information is readily available renders the analysis much more straightforward in the H→ZZ→4e¡° channels. | |||||||||||||||||||||||||||||||||||||||||||||
> > | The only channel where all information about the boson is available in the detector is the H→ZZ→4e“ where the leptons are electrons or muons (and not neutrinos). This does not mean that other channels cannot contribute to this picture (like H→WW→e“νe“ν or H→γγ), but the fact that all information is readily available renders the analysis much more straightforward in the H→ZZ→4e“ channels. | |||||||||||||||||||||||||||||||||||||||||||||
The following is an illustration from http://arxiv.org/pdf/1208.4018.pdf![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Since the spin-parity tensor structure is much more complicated than the scalar structure of the couplings, much more data will be needed to perform precise measurements. For the moment the strategy is to perform hypotheses tests between the 0+ (SMH-like) hypothesis and other, alternative, hypotheses such as pseudoscalar (0-) bosons (or admixtures), spin-2 bosons, like the graviton, (2+), or exotic spin-1 bosons. As an example, the following figure illustrates the differences in one observable of the ZZ (left), WW (middle), and γγ (right) decays for different JCP hypotheses (0+ red circles, 0- magenta squares, and two types of spin-2 particle in green and blue): | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Intro 2: Datacards | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 240 to 239 | ||||||||||||||||||||||||||||||||||||||||||||||
CMS_vbfhinv_zvv_stat gmN 12 - 8.4899 - - - - - | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | The gamma distribution, in the form we use it, is characterized by two parameters: N, the number of events in the control region (or in MC), and α the weight by which each event is multiplied to determine the expected yield in the signal region (if the weights are not all the same, one can take α to be the average weight). By definition, the prediciton in the signal region is | |||||||||||||||||||||||||||||||||||||||||||||
> > | The gamma distribution, in the form we use it, is characterized by two parameters: N, the number of events in the control region (or in MC), and α the weight by which each event is multiplied to determine the expected yield in the signal region (if the weights are not all the same, one can take α to be the average weight). By definition, the prediciton in the signal region is | |||||||||||||||||||||||||||||||||||||||||||||
In the datacards, the uncertainty is encoded writing gmN immediately followed by the number of events N, and then putting the factor α in the column corresponding to the process. So, in this case N=12, α=8.5 and the expected yield is about 102 events. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 270 to 269 | ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | So, keeping only the relevant decimals, the total expected background is 340 ¡À 18(stat) ¡À 61 (syst). | |||||||||||||||||||||||||||||||||||||||||||||
> > | So, keeping only the relevant decimals, the total expected background is 340 ± 18(stat) ± 61 (syst). | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> A shape analysis using templates: H→ττ | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 342 to 344 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% /// root -l root -l htt_mt.input_8TeV.root /// Attaching file htt_mt.input_8TeV.root as _file0... | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | _file0->cd("muTau_vbf"); QCD->Draw("HIST"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | _file0->cd("muTau_vbf"); QCD->Draw("HIST"); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
We can also check the normalization of the histogram, to see that it matches the expected yield in the datacard as it should. %CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | cout << QCD->Integral() << endl; /// ----> 39.2563 | |||||||||||||||||||||||||||||||||||||||||||||
> > | cout << QCD->Integral() << endl; /// ----> 39.2563 | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% Note that underflow and overflow bins are not used by combine. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 359 to 361 | ||||||||||||||||||||||||||||||||||||||||||||||
Similarly we can look at how the 199 events observed in data are distributed: the first rule again will give muTau_vbf/data_obs
%CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | cout << data_obs->Integral() << endl; /// ----> 199 data_obs->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | cout << data_obs->Integral() << endl; /// ----> 199 data_obs->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
We can also at how the signal looks like for different Higgs boson mass hypotheses; we will use the useful Form("pattern", args) method of ROOT to create the proper string with the mass written inside.
%CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | qqH110->Draw("HIST"); for (int m = 110; m <= 140; m += 10) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); h->SetLineColor(206+4*(m-110)/10); h->SetLineWidth(2); h->Draw("HIST SAME"); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | qqH110->Draw("HIST"); for (int m = 110; m <= 140; m += 10) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); h->SetLineColor(206+4*(m-110)/10); h->SetLineWidth(2); h->Draw("HIST SAME"); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% for (int m = 110; m <= 140; m += 5) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); printf("for mH=%3d, signal mean = %5.1f, rms = %5.1f\n", m, h->GetMean(), h->GetRMS()); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% for (int m = 110; m <= 140; m += 5) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); printf("for mH=%3d, signal mean = %5.1f, rms = %5.1f\n", m, h->GetMean(), h->GetRMS()); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
for mH=110, signal mean = 100.5, rms = 15.1 for mH=115, signal mean = 104.6, rms = 16.8 for mH=120, signal mean = 108.8, rms = 17.5 | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 402 to 406 | ||||||||||||||||||||||||||||||||||||||||||||||
This tells us that there is a shape uncertainty called CMS_scale_t_mutau_8TeV which affects the signals and the ZTT background (this background is determined embedding simulated taus into Z→μμ events, so any data/mc discrepancy on the energy response for taus will affect the prediction for it just like it does for the simulated taus in signal MC) | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | For shape analyses using templates, the effect of a shape-affecting systematical uncertainty is encoded in the datacard by providing two additional templates corresponding to ¡ÀNσ shifts of the corresponding parameter. For each process affected by this uncertainty, the column in the datacard will contain a value 1/N. In this case, as often, ¡À1σ shapes have been used, and so the columns of the datacard contain a value of 1. | |||||||||||||||||||||||||||||||||||||||||||||
> > | For shape analyses using templates, the effect of a shape-affecting systematical uncertainty is encoded in the datacard by providing two additional templates corresponding to ±Nσ shifts of the corresponding parameter. For each process affected by this uncertainty, the column in the datacard will contain a value 1/N. In this case, as often, ±1σ shapes have been used, and so the columns of the datacard contain a value of 1. | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | The associated shapes are determined using the last column of the shapes line, which for this datacard reads $CHANNEL/$PROCESS_$SYSTEMATIC , where $SYSTEMATIC is taken to be the name of the uncertainty in the card, and the postixes Up and Down are added to get the up and down templates. Let's see how much does the mean of the expected ZTT distribution change for the ¡À1σ shifts, and also look at the shapes in ROOT. We can also see that the tau energy scale does not affect only the shape of the background, but also the normalization (probably because if taus are reconstructed as more energetic, they have a higher probability of passing the pT cuts or other kinematic cuts of the analsysis) | |||||||||||||||||||||||||||||||||||||||||||||
> > | The associated shapes are determined using the last column of the shapes line, which for this datacard reads $CHANNEL/$PROCESS_$SYSTEMATIC , where $SYSTEMATIC is taken to be the name of the uncertainty in the card, and the postixes Up and Down are added to get the up and down templates. Let's see how much does the mean of the expected ZTT distribution change for the ±1σ shifts, and also look at the shapes in ROOT. We can also see that the tau energy scale does not affect only the shape of the background, but also the normalization (probably because if taus are reconstructed as more energetic, they have a higher probability of passing the pT cuts or other kinematic cuts of the analsysis) | |||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% /// root -l htt_mt.input_8TeV.root /// Attaching file htt_mt.input_8TeV.root as _file0... | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | _file0->cd("muTau_vbf"); TH1F nominal = (TH1F) gDirectory->Get("ZTT"); TH1F tau_up = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVUp"); TH1F tau_dn = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVDown"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | _file0->cd("muTau_vbf"); TH1F nominal = (TH1F) gDirectory->Get("ZTT"); TH1F tau_up = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVUp"); TH1F tau_dn = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVDown"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | printf("Nominal shape: mean = %6.2f (norm: %7.3f events)\n", nominal->GetMean(), nominal->Integral()); // --> Nominal shape: mean = 88.03 (norm: 94.286 events) printf("Up shape: mean = %6.2f (norm: %7.3f events)\n", tau_up->GetMean(), tau_up->Integral()); // --> Up shape: mean = 88.86 (norm: 97.172 events) printf("Down shape: mean = %6.2f (norm: %7.3f events)\n", tau_dn->GetMean(), tau_dn->Integral()); // --> Down shape: mean = 87.20 (norm: 91.995 events) | |||||||||||||||||||||||||||||||||||||||||||||
> > | printf("Nominal shape: mean = %6.2f (norm: %7.3f events)\n", nominal->GetMean(), nominal->Integral()); // --> Nominal shape: mean = 88.03 (norm: 94.286 events) printf("Up shape: mean = %6.2f (norm: %7.3f events)\n", tau_up->GetMean(), tau_up->Integral()); // --> Up shape: mean = 88.86 (norm: 97.172 events) printf("Down shape: mean = %6.2f (norm: %7.3f events)\n", tau_dn->GetMean(), tau_dn->Integral()); // --> Down shape: mean = 87.20 (norm: 91.995 events) | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | nominal->SetLineWidth(2); nominal->SetLineColor(kBlack); tau_up->SetLineWidth(2); tau_up->SetLineColor(kBlue); tau_dn->SetLineWidth(2); tau_dn->SetLineColor(kRed); nominal->Draw("HIST"); tau_up->Draw("HIST SAME"); tau_dn->Draw("HIST SAME"); nominal->GetYaxis()->SetRangeUser(0,45); | |||||||||||||||||||||||||||||||||||||||||||||
> > | nominal->SetLineWidth(2); nominal->SetLineColor(kBlack); tau_up->SetLineWidth(2); tau_up->SetLineColor(kBlue); tau_dn->SetLineWidth(2); tau_dn->SetLineColor(kRed); nominal->Draw("HIST"); tau_up->Draw("HIST SAME"); tau_dn->Draw("HIST SAME"); nominal->GetYaxis()->SetRangeUser(0,45); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
A parametric shape analysis: H→γγ | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 498 to 504 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% TFile *fDat = TFile::Open("hgg.inputbkgdata_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooAbsData *data = cms_hgg_workspace->data("roohist_data_mass_cat0"); data->Print("") // --> RooDataHist::roohist_data_mass_cat0[CMS_hgg_mass] = 160 bins (1449 weights) | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooAbsData *data = cms_hgg_workspace->data("roohist_data_mass_cat0"); data->Print("") // --> RooDataHist::roohist_data_mass_cat0[CMS_hgg_mass] = 160 bins (1449 weights) | |||||||||||||||||||||||||||||||||||||||||||||
// so, we have a binned dataset, whose variable is called CMS_hgg_mass: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooRealVar *mass = cms_hgg_workspace->var("CMS_hgg_mass"); mass->Print(""); | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooRealVar *mass = cms_hgg_workspace->var("CMS_hgg_mass"); mass->Print(""); | |||||||||||||||||||||||||||||||||||||||||||||
// RooRealVar::CMS_hgg_mass = 140 L(100 - 180) // we can make a plot of the dataset with the following | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooPlot *plot = mass->frame(); data->plotOn(plot); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooPlot *plot = mass->frame(); data->plotOn(plot); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 518 to 524 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); ggH->Print(""); // --> RooAddPdf::hggpdfrel_ggh_cat0[ hist_func_frac_g0_ggh_cat0 * hgg_gaus_g0_ggh_cat0 + hist_func_frac_g1_ggh_cat0 * hgg_gaus_g1_ggh_cat0 + const_func_frac_g2_ggh_cat0 * hgg_gaus_g2_ggh_cat0 ] = 0.604097 // this appears to be a linear combination of multiple gaussian pdfs (hgg_gaus_g0_ggh_cat0, hgg_gaus_g1_ggh_cat0, ...) with coefficents hist_func_frac_g0_ggh_cat0, hist_func_frac_g1_ggh_cat0 | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); ggH->Print(""); // --> RooAddPdf::hggpdfrel_ggh_cat0[ hist_func_frac_g0_ggh_cat0 * hgg_gaus_g0_ggh_cat0 + hist_func_frac_g1_ggh_cat0 * hgg_gaus_g1_ggh_cat0 + const_func_frac_g2_ggh_cat0 * hgg_gaus_g2_ggh_cat0 ] = 0.604097 // this appears to be a linear combination of multiple gaussian pdfs (hgg_gaus_g0_ggh_cat0, hgg_gaus_g1_ggh_cat0, ...) with coefficents hist_func_frac_g0_ggh_cat0, hist_func_frac_g1_ggh_cat0 | |||||||||||||||||||||||||||||||||||||||||||||
// let's get the list of the parameters that describe it | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooArgSet *params = ggH->getParameters(*data); params->Print(""); // --> RooArgSet::parameters = (CMS_hgg_globalscale,CMS_hgg_nuissancedeltamcat0,CMS_hgg_nuissancedeltasmearcat0,MH) | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooArgSet *params = ggH->getParameters(*data); params->Print(""); // --> RooArgSet::parameters = (CMS_hgg_globalscale,CMS_hgg_nuissancedeltamcat0,CMS_hgg_nuissancedeltasmearcat0,MH) | |||||||||||||||||||||||||||||||||||||||||||||
// MH is a special parameter, which combine and text2workspace set to the Higgs mass hypothesis. // now since we're just looking into the input workspace, we can set it by hand | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("MH")->setVal(125.7); | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("MH")->setVal(125.7); | |||||||||||||||||||||||||||||||||||||||||||||
// Now we can make a plot of the pdf, and show also the contributions from the different gaussians from which it is composed | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame(); ggH->plotOn(plot); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g0_ggh_cat0"), RooFit::LineColor(kRed)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g1_ggh_cat0"), RooFit::LineColor(209)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g2_ggh_cat0"), RooFit::LineColor(222)); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame(); ggH->plotOn(plot); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g0_ggh_cat0"), RooFit::LineColor(kRed)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g1_ggh_cat0"), RooFit::LineColor(209)); ggH->plotOn(plot, RooFit::Components("hgg_gaus_g2_ggh_cat0"), RooFit::LineColor(222)); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 555 to 561 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); RooAbsReal *ggH_norm = wsig_8TeV->function("hggpdfrel_ggh_cat0_norm"); cout << ggH_norm->getVal()*19620.0000 << endl; // --> 12.3902 | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); RooAbsReal *ggH_norm = wsig_8TeV->function("hggpdfrel_ggh_cat0_norm"); cout << ggH_norm->getVal()*19620.0000 << endl; // --> 12.3902 | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
This approach can also be used to make a background freely floating, by associating to them as normalization term a floating RooRealVar. However, this can be achieved in a more transparent way by putting instead a normalization uncertainty on that background using a flat pdf: e.g. to leave the background floating between 50% and 200% of its input prediction, a lnU systematic can be used with κ=2 ( lnU has a syntax like lnN , but produces a uniform pdf between 1/κ and κ rather than a log-normal; see SWGuideHiggsAnalysisCombinedLimit![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 569 to 575 | ||||||||||||||||||||||||||||||||||||||||||||||
The part of the H→γγ datacard related to the systematics starts with many lines of log-normals that should already be familiar, except possibly for the notation with two numbers separated by a slash (e.g. 0.950/1.099). This notation is used for asymmetrical uncertainties: a log-normal with 0.950/1.099 means that at -1σ the yield is scaled down by a factor 0.95, while at +1σ the yield is scaled up by a factor 1.099. The last part of the datacard contains some lines that use a different syntax, e.g. | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
CMS_hgg_globalscale param 0.0 0.004717 | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | These lines encode uncertainties on the parameters of the signal and background pdfs. The example line quoted here informs text2workspace that the parameter CMS_hgg_globalscale is to be assigned a Gaussian uncertainty of ¡À0.004717 around its mean value of zero (0.0). | |||||||||||||||||||||||||||||||||||||||||||||
> > | These lines encode uncertainties on the parameters of the signal and background pdfs. The example line quoted here informs text2workspace that the parameter CMS_hgg_globalscale is to be assigned a Gaussian uncertainty of ±0.004717 around its mean value of zero (0.0). | |||||||||||||||||||||||||||||||||||||||||||||
The effect can be visualized from RooFit, e.g. %CODE{"cpp"}% TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("MH")->setVal(125.7); RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0"); | |||||||||||||||||||||||||||||||||||||||||||||
// prepare the canvas | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame(); | |||||||||||||||||||||||||||||||||||||||||||||
// plot nominal pdf | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ggH->plotOn(plot, RooFit::LineColor(kBlack)); | |||||||||||||||||||||||||||||||||||||||||||||
> > | ggH->plotOn(plot, RooFit::LineColor(kBlack)); | |||||||||||||||||||||||||||||||||||||||||||||
// plot minus 3 sigma pdf | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("CMS_hgg_globalscale")->setVal(-3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kBlue)); | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("CMS_hgg_globalscale")->setVal(-3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kBlue)); | |||||||||||||||||||||||||||||||||||||||||||||
// plot plus 3 sigma pdf | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("CMS_hgg_globalscale")->setVal(+3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kRed)); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("CMS_hgg_globalscale")->setVal(+3*0.004717); ggH->plotOn(plot, RooFit::LineColor(kRed)); plot->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 612 to 619 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% TFile *fDat = TFile::Open("hgg.inputbkgdata_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | cout << cms_hgg_workspace->data("roohist_data_mass_cat0")->sumEntries() << endl; | |||||||||||||||||||||||||||||||||||||||||||||
> > | cout << cms_hgg_workspace->data("roohist_data_mass_cat0")->sumEntries() << endl; | |||||||||||||||||||||||||||||||||||||||||||||
// .... and so on %ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 636 to 642 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% TFile *fSig = TFile::Open("hgg.inputsig_8TeV_MVA.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | wsig_8TeV->var("MH")->setVal(125.7); cout << wsig_8TeV->function("hggpdfrel_ggh_cat0_norm")->getVal()*19620.0000 << endl; cout << wsig_8TeV->function("hggpdfrel_vbf_cat0_norm")->getVal()*19620.0000 + wsig_8TeV->function("hggpdfrel_wh_cat0_norm")->getVal()*12753.0000 + wsig_8TeV->function("hggpdfrel_zh_cat0_norm")->getVal()*6867.0000 << endl; // .... and so on | |||||||||||||||||||||||||||||||||||||||||||||
> > | wsig_8TeV->var("MH")->setVal(125.7); cout << wsig_8TeV->function("hggpdfrel_ggh_cat0_norm")->getVal()*19620.0000 << endl; cout << wsig_8TeV->function("hggpdfrel_vbf_cat0_norm")->getVal()*19620.0000 + wsig_8TeV->function("hggpdfrel_wh_cat0_norm")->getVal()*12753.0000 + wsig_8TeV->function("hggpdfrel_zh_cat0_norm")->getVal()*6867.0000 << endl; // .... and so on | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
| ||||||||||||||||||||||||||||||||||||||||||||||
Line: 703 to 709 | ||||||||||||||||||||||||||||||||||||||||||||||
The interesting difference however is the observed dataset: %CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | gSystem->Load("libHiggsAnalysisCombinedLimit.so") | |||||||||||||||||||||||||||||||||||||||||||||
> > | gSystem->Load("libHiggsAnalysisCombinedLimit.so") | |||||||||||||||||||||||||||||||||||||||||||||
gFile = TFile::Open("hzz4l_4muS_8TeV_0.input.root"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooAbsData *data = w->data("data_obs"); data->Print("") // --> RooDataSet::data_obs[CMS_zz4l_mass,melaLD,CMS_zz4l_PToverM] = 13 entries | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooAbsData *data = w->data("data_obs"); data->Print("") // --> RooDataSet::data_obs[CMS_zz4l_mass,melaLD,CMS_zz4l_PToverM] = 13 entries | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% In the Hγγ the dataset was a RooDataHist with one entry per bin, each entry having a weight equal to the number of events in that bin. Here instead it's a RooDataSet, with one entry per event. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 715 to 721 | ||||||||||||||||||||||||||||||||||||||||||||||
Let's explore its contents, both building some binned 1D projections and by making a 2D scatter plot. %CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | w->var("CMS_zz4l_mass")->Print(""); // --> RooRealVar::CMS_zz4l_mass = 114.157 L(104 - 142) B(200) RooPlot *mass = w->var("CMS_zz4l_mass")->frame(RooFit::Bins(19)); data->plotOn(mass); mass->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | w->var("CMS_zz4l_mass")->Print(""); // --> RooRealVar::CMS_zz4l_mass = 114.157 L(104 - 142) B(200) RooPlot *mass = w->var("CMS_zz4l_mass")->frame(RooFit::Bins(19)); data->plotOn(mass); mass->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | RooPlot *mela = w->var("melaLD")->frame(RooFit::Bins(10)); data->plotOn(mass); mass->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | RooPlot *mela = w->var("melaLD")->frame(RooFit::Bins(10)); data->plotOn(mass); mass->Draw(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TGraph points(data->numEntries()); for (int i = 0; i < data->numEntries(); ++i) { const RooArgSet *ev = data->get(i); points->SetPoint(i, ev->getRealValue("CMS_zz4l_mass"), ev->getRealValue("melaLD")); } points->SetMarkerStyle(20); points->GetXaxis()->SetTitle("m_{4l} (GeV)") points->GetYaxis()->SetTitle("K_{D}") .x tdrstyle.cc points->Draw("AP"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TGraph points(data->numEntries()); for (int i = 0; i < data->numEntries(); ++i) { const RooArgSet *ev = data->get(i); points->SetPoint(i, ev->getRealValue("CMS_zz4l_mass"), ev->getRealValue("melaLD")); } points->SetMarkerStyle(20); points->GetXaxis()->SetTitle("m_{4l} (GeV)") points->GetYaxis()->SetTitle("K_{D}") .x tdrstyle.cc points->Draw("AP"); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
Combining datacards | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 846 to 858 | ||||||||||||||||||||||||||||||||||||||||||||||
Now we can open the hww2j.root file in ROOT and see its contents. Start a root prompt with the command root.exe -b -l -n and | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // Load the libraries needed to interpret the likelihood gSystem->Load("libHiggsAnalysisCombinedLimit.so"); // Open the file gFile = TFile::Open("hww2j.root"); // Retrieve the workspace from the file RooWorkspace *w = (RooWorkspace *) gFile->Get("w"); // Retrieve the ModelConfig from the workspace RooStats::ModelConfig *mc = (RooStats::ModelConfig *) w->genobj("ModelConfig"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // Load the libraries needed to interpret the likelihood gSystem->Load("libHiggsAnalysisCombinedLimit.so"); // Open the file gFile = TFile::Open("hww2j.root"); // Retrieve the workspace from the file RooWorkspace *w = (RooWorkspace *) gFile->Get("w"); // Retrieve the ModelConfig from the workspace RooStats::ModelConfig *mc = (RooStats::ModelConfig *) w->genobj("ModelConfig"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // Now we get the list of parameters of the workspace mc->GetParametersOfInterest()->Print("V"); // the output will be // 1) RooRealVar:: r = 1 | |||||||||||||||||||||||||||||||||||||||||||||
> > | // Now we get the list of parameters of the workspace mc->GetParametersOfInterest()->Print("V"); // the output will be // 1) RooRealVar:: r = 1 | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // Now we get the expected signal yield for ggH and DYTT cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // the output will be 1.44, and 1.879, as expected | |||||||||||||||||||||||||||||||||||||||||||||
> > | // Now we get the expected signal yield for ggH and DYTT cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // the output will be 1.44, and 1.879, as expected | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // Now we change the value of r to be 2.0... we'd expect the ggH yield to grow to 2.88 and the DYTT to remain the same w->var("r")->setVal(2.0); cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; // ----> 2.88 cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // ----> 1.879 | |||||||||||||||||||||||||||||||||||||||||||||
> > | // Now we change the value of r to be 2.0... we'd expect the ggH yield to grow to 2.88 and the DYTT to remain the same w->var("r")->setVal(2.0); cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; // ----> 2.88 cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // ----> 1.879 | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // We can also get the number of observed events, that should match the "observation" line in the datacard cout << w->data("data_obs")->sumEntries() << endl; // ----> 24 | |||||||||||||||||||||||||||||||||||||||||||||
> > | // We can also get the number of observed events, that should match the "observation" line in the datacard cout << w->data("data_obs")->sumEntries() << endl; // ----> 24 | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
One-dimensional fits and profile likelihood | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 880 to 894 | ||||||||||||||||||||||||||||||||||||||||||||||
The output file from combine will contain a tree with a single entry, corresponding to this value. | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // root.exe -b -l higgsCombineHWW_FIT.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_FIT.MultiDimFit.mH125.7.root as _file0... root [1] limit->Scan("r"); //************************ //* Row * r * //************************ //* 0 * 0.5986483 * //************************ | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // root.exe -b -l higgsCombineHWW_FIT.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_FIT.MultiDimFit.mH125.7.root as _file0... root [1] limit->Scan("r"); //************************ //* Row * r * //************************ //* 0 * 0.5986483 * //************************ | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
Now, to compare different values of μ, we will use the negative log of the profile likelihood function λ | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 898 to 913 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% // root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | limit->GetEntries() // (const Long64_t)101 | |||||||||||||||||||||||||||||||||||||||||||||
> > | limit->GetEntries() // (const Long64_t)101 | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | limit->Scan("r:2*deltaNLL") //************************************ //* Row * r * 2*deltaNL * //************************************ //* 0 * 0.5986473 * 0 * //* 1 * 0.0149999 * 5.9730005 * //* 2 * 0.0450000 * 5.3444776 * //* 3 * 0.0750000 * 4.7531828 * //* 4 * 0.1049999 * 4.1995768 * //* 5 * 0.1350000 * 3.6845009 * | |||||||||||||||||||||||||||||||||||||||||||||
> > | limit->Scan("r:2*deltaNLL") //************************************ //* Row * r * 2*deltaNL * //************************************ //* 0 * 0.5986473 * 0 * //* 1 * 0.0149999 * 5.9730005 * //* 2 * 0.0450000 * 5.3444776 * //* 3 * 0.0750000 * 4.7531828 * //* 4 * 0.1049999 * 4.1995768 * //* 5 * 0.1350000 * 3.6845009 * | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | limit->Draw("2*deltaNLL:r"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | limit->Draw("2*deltaNLL:r"); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
The quality of the default plot looks awful, so we will retrieve that information as a TGraph object and make it a bit more pretty; also, we will zoom into the region closer to the minimum which is more interesting. | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // Draw limit->Draw("2*deltaNLL:r"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // Draw limit->Draw("2*deltaNLL:r"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // Make a copy of the graph TGraph *graph = gROOT->FindObject("Graph")->Clone(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // Make a copy of the graph TGraph *graph = gROOT->FindObject("Graph")->Clone(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // sort the points for increasing values of r graph->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // sort the points for increasing values of r graph->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // set the style graph->SetLineWidth(2); graph->SetLineColor(214); graph->SetMarkerStyle(7); graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8); c1->SetGridy(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // set the style graph->SetLineWidth(2); graph->SetLineColor(214); graph->SetMarkerStyle(7); graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8); c1->SetGridy(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // draw as a smooth line plus the points graph->Draw("APC"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // draw as a smooth line plus the points graph->Draw("APC"); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Deleted: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Now, from the scan you could read the boundaries of the ¡À1σ interval by looking at where q(μ) crosses the horizontal line at 1. However, since that's a very common task, for models with just one parameter combine can do that for you using the MaxLikelihoodFit module | |||||||||||||||||||||||||||||||||||||||||||||
combine -M MaxLikelihoodFit comb_hww.root -m 125.7 -n HWW_FIT_MU --justFit --robustFit=1 | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1027 to 1045 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp" }% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% /// root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root as _file0... // Attaching file higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graph = gROOT->FindObject("Graph")->Clone(); graph->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% /// root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root as _file0... // Attaching file higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graph = gROOT->FindObject("Graph")->Clone(); graph->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | _file1->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graphA = gROOT->FindObject("Graph")->Clone(); graphA->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | _file1->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graphA = gROOT->FindObject("Graph")->Clone(); graphA->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | graph->SetLineWidth(3); graph->SetLineColor(1); graphA->SetLineWidth(1); graphA->SetLineColor(4); graphA->SetLineStyle(2); graph->Draw("AC"); graphA->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | graph->SetLineWidth(3); graph->SetLineColor(1); graphA->SetLineWidth(1); graphA->SetLineColor(4); graphA->SetLineStyle(2); graph->Draw("AC"); graphA->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8); | |||||||||||||||||||||||||||||||||||||||||||||
> > | graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1101 to 1121 | ||||||||||||||||||||||||||||||||||||||||||||||
A simple plot from the output of the bidimensional scan can be done with %CODE{"cpp"}% // plot the likelihood in the 20x20 grid, selecting the interesting region of small deltaNLL | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // make the plot less horrible gStyle->SetOptStat(0); h2->SetContour(200); h2->SetTitle("Profile Likelihood in 2D;#mu_{ggH,ttH};#mu_{VBF,VH}"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // make the plot less horrible gStyle->SetOptStat(0); h2->SetContour(200); h2->SetTitle("Profile Likelihood in 2D;#mu_{ggH,ttH};#mu_{VBF,VH}"); | |||||||||||||||||||||||||||||||||||||||||||||
/// Now we add the best fit point | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | limit->Draw("RV:RF","deltaNLL == 0","P SAME"); TGraph *bestFit = gROOT->FindObject("Graph"); bestFit->SetMarkerStyle(34); bestFit->SetMarkerSize(3); c1->Update(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | limit->Draw("RV:RF","deltaNLL == 0","P SAME"); TGraph *bestFit = gROOT->FindObject("Graph"); bestFit->SetMarkerStyle(34); bestFit->SetMarkerSize(3); c1->Update(); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Contours for different confidence levels can be constructed starting from the 2d likelihood scan, using the thresholds for a chisquare with two degrees of freedom (pdf review![]() ROOT::Math::chisquared_quantile(cl,ndf) can be used to compute the thresholds for each CL. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1133 to 1153 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root contours2D.cxx // Attaching file higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root as _file0... // Processing contours2D.cxx... | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | // before we had done // limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | // before we had done // limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); | |||||||||||||||||||||||||||||||||||||||||||||
// now we do contour2D("RF",20,-1.,3., "RV",20,-2.,5., 1.,1.); %ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1149 to 1170 | ||||||||||||||||||||||||||||||||||||||||||||||
combine -M MultiDimFit rVrF_hgg.root -m 125.7 --algo=grid --points=50 --setPhysicsModelParameterRanges RF=-1,3 -P RF --floatOtherPOI=1 -n HGG_RVRF_SCAN_1D_RF_PROF | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_1D_RF_* /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_FIX.MultiDimFit.mH125.7.root as _file0... /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_PROF.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphFix = gROOT->FindObject("Graph")->Clone(); graphFix->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_1D_RF_* /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_FIX.MultiDimFit.mH125.7.root as _file0... /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_PROF.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphFix = gROOT->FindObject("Graph")->Clone(); graphFix->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | _file1->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphProf = gROOT->FindObject("Graph")->Clone(); graphProf->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | _file1->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphProf = gROOT->FindObject("Graph")->Clone(); graphProf->Sort(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | graphProf->SetLineWidth(3); graphProf->Draw("AC"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | graphProf->SetLineWidth(3); graphProf->Draw("AC"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | graphFix->SetLineWidth(2); graphFix->SetLineStyle(9); graphFix->SetLineColor(206); graphFix->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | graphFix->SetLineWidth(2); graphFix->SetLineStyle(9); graphFix->SetLineColor(206); graphFix->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1219 to 1242 | ||||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
A macro to make all the plots including the summary one is below: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotall_rvrf.cxx -q void plotall_rvrf() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("rvrf_all.root", "RECREATE"); const int decays = 5; const char *decay[decays] = { "hww", "hzz", "hgg", "hbb", "htt" }; const int color[decays] = { 4 , 2 , 209 , 67 , 221 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_RVRF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("RF",60,-1,3,"RV",60,-2,5,1.0,1.0,fOut,Form("rvrf_%s",decay[i])); c1->Print(Form("output_rvrf_%s.png",decay[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#mu_{ggH,ttH};#mu_{VBF,VH}",60,-1.,3.,60,-2.,5.); frame->Draw("AXIS"); for (int i = 0; i < decays; ++i) { TGraph *best = fOut->Get(Form("rvrf_%s_best",decay[i])); best->SetMarkerColor(color[i]); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); TList *c68 = fOut->Get(Form("rvrf_%s_c68",decay[i])); styleMultiGraph(c68, /*color=*/color[i], /*width=*/3, /*style=*/1); c68->Draw("L SAME"); } TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_rvrf_all.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotall_rvrf.cxx -q void plotall_rvrf() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("rvrf_all.root", "RECREATE"); const int decays = 5; const char *decay[decays] = { "hww", "hzz", "hgg", "hbb", "htt" }; const int color[decays] = { 4 , 2 , 209 , 67 , 221 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_RVRF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("RF",60,-1,3,"RV",60,-2,5,1.0,1.0,fOut,Form("rvrf_%s",decay[i])); c1->Print(Form("output_rvrf_%s.png",decay[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#mu_{ggH,ttH};#mu_{VBF,VH}",60,-1.,3.,60,-2.,5.); frame->Draw("AXIS"); for (int i = 0; i < decays; ++i) { TGraph *best = fOut->Get(Form("rvrf_%s_best",decay[i])); best->SetMarkerColor(color[i]); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); TList *c68 = fOut->Get(Form("rvrf_%s_c68",decay[i])); styleMultiGraph(c68, /*color=*/color[i], /*width=*/3, /*style=*/1); c68->Draw("L SAME"); } TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_rvrf_all.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1244 to 1266 | ||||||||||||||||||||||||||||||||||||||||||||||
The loop-induced effective couplings of the Higgs boson to gluons and to photons are computed in terms of κV, κf by assuming that only the SM particles participate to the loop.
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Cross sections and partial all widths scale quadratically with the couplings, so in this model we'll have e.g. | |||||||||||||||||||||||||||||||||||||||||||||
> > | Cross sections and partial all widths scale quadratically with the couplings, so in this model we'll have e.g. | |||||||||||||||||||||||||||||||||||||||||||||
In this model we assume that there are no BSM decays, so where we have neglected the γγ and Zγ decays that contribute only to the few per-mille level. We can then dividing by the SM width, so that we can express a scale factor for the total width κH in terms of SM BRs Since in the SM for a 125.7 GeV Higgs we have BR(H→VV) = 0.254, to a good approximation Exercises for the reader | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
%TWISTY{ mode="div" showlink="Answer (don't look until you've tried it!)" | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1261 to 1283 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp"
}%
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
It is instructive to note that in this benchmark model the search for VH, H→bb provides exactly the same kind of measurement as ggH, H→ZZ (or ggH, H→VV). The same is true also for another channel which is often intuitively associated with fermion couplings, VBF H→ττ.
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1280 to 1302 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp"
}%
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Among the channels considered, the ones with the strongest dependency on κV2 and κf2 are respectively VBF H→γγ and ttH, H→bb. More in general, you can deduce that the strongest dependency is achieved when considering final states with the same coupling both at production and decay (VBF and VH production with H→γγ,VV; ggH and ttH production with H→bb, τtau;), while the channels with mixed productions and decays have smaller coefficients since there is a partial cancellation between the κf2 at numerator and the ¾κf2 at denominator (from κH2). | |||||||||||||||||||||||||||||||||||||||||||||
> > | Among the channels considered, the ones with the strongest dependency on κV2 and κf2 are respectively VBF H→γγ and ttH, H→bb. More in general, you can deduce that the strongest dependency is achieved when considering final states with the same coupling both at production and decay (VBF and VH production with H→γγ,VV; ggH and ttH production with H→bb, τtau;), while the channels with mixed productions and decays have smaller coefficients since there is a partial cancellation between the κf2 at numerator and the ¾κf2 at denominator (from κH2). | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> Implementation in combine | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1351 to 1375 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp" }% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_hww.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_hww() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_hww.root", "RECREATE"); const int states = 3; const char *state[states] = { "hww", "hww0j", "hww2j" }; for (int i = 0; i < states; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D_ASIMOV.MultiDimFit.mH125.7.root", state[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",state[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_hww.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_hww() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_hww.root", "RECREATE"); const int states = 3; const char *state[states] = { "hww", "hww0j", "hww2j" }; for (int i = 0; i < states; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D_ASIMOV.MultiDimFit.mH125.7.root", state[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",state[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TList c68_comb = (TList) fOut->Get("kvkf_hww_c68"); styleMultiGraph(c68_comb, /*color=*/1, /*width=*/3, /*style=*/1); fillMultiGraph(c68_comb, /*color=*/kGray, /*style=*/1001); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TList c68_comb = (TList) fOut->Get("kvkf_hww_c68"); styleMultiGraph(c68_comb, /*color=*/1, /*width=*/3, /*style=*/1); fillMultiGraph(c68_comb, /*color=*/kGray, /*style=*/1001); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TList c68_0j = (TList) fOut->Get("kvkf_hww0j_c68"); styleMultiGraph(c68_0j, /*color=*/4, /*width=*/3, /*style=*/1); fillMultiGraph(c68_0j, /*color=*/4, /*style=*/3005); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TList c68_0j = (TList) fOut->Get("kvkf_hww0j_c68"); styleMultiGraph(c68_0j, /*color=*/4, /*width=*/3, /*style=*/1); fillMultiGraph(c68_0j, /*color=*/4, /*style=*/3005); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TList c68_2j = (TList) fOut->Get("kvkf_hww2j_c68"); fillMultiGraph(c68_2j, /*color=*/2, /*style=*/3004); styleMultiGraph(c68_2j, /*color=*/2, /*width=*/3, /*style=*/1); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TList c68_2j = (TList) fOut->Get("kvkf_hww2j_c68"); fillMultiGraph(c68_2j, /*color=*/2, /*style=*/3004); styleMultiGraph(c68_2j, /*color=*/2, /*width=*/3, /*style=*/1); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | c68_2j->Draw("F SAME"); c68_0j->Draw("F SAME"); c68_comb->Draw("F SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | c68_2j->Draw("F SAME"); c68_0j->Draw("F SAME"); c68_comb->Draw("F SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | c68_2j->Draw("L SAME"); c68_0j->Draw("L SAME"); c68_comb->Draw("L SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | c68_2j->Draw("L SAME"); c68_0j->Draw("L SAME"); c68_comb->Draw("L SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TGraph *best = fOut->Get("kvkf_hww_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TGraph *best = fOut->Get("kvkf_hww_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); TLegend *leg = new TLegend(0.20,0.80,0.45,0.95); leg->SetTextFont(42); leg->SetShadowColor(0); leg->SetFillColor(0); leg->AddEntry(best, "Comb.", "P"); leg->AddEntry(c68_comb->At(0), "(68% CL)", "LF"); leg->AddEntry(c68_0j->At(0), "0 jets", "LF"); leg->AddEntry(c68_2j->At(0), "2 jets", "LF"); leg->Draw(); c1->Print("output_kvkf_hww.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); TLegend *leg = new TLegend(0.20,0.80,0.45,0.95); leg->SetTextFont(42); leg->SetShadowColor(0); leg->SetFillColor(0); leg->AddEntry(best, "Comb.", "P"); leg->AddEntry(c68_comb->At(0), "(68% CL)", "LF"); leg->AddEntry(c68_0j->At(0), "0 jets", "LF"); leg->AddEntry(c68_2j->At(0), "2 jets", "LF"); leg->Draw(); c1->Print("output_kvkf_hww.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1401 to 1425 | ||||||||||||||||||||||||||||||||||||||||||||||
showimgleft="/pub/TWiki/TWikiDocGraphics/toggleopen-small.gif" hideimgleft="/pub/TWiki/TWikiDocGraphics/toggleclose-small.gif" class="twikiHelp" | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | }% Γtot/ΓtotSM = 0.566¡¤κb2 + 0.2541¡¤κV2 + 0.0851¡¤κg2 + 0.0623¡¤κτ2 + 0.0285¡¤κt2 + 0.00388¡¤κγ2 | |||||||||||||||||||||||||||||||||||||||||||||
> > | }% Γtot/ΓtotSM = 0.566·κb2 + 0.2541·κV2 + 0.0851·κg2 + 0.0623·κτ2 + 0.0285·κt2 + 0.00388·κγ2 | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> Implementation in combine | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1426 to 1450 | ||||||||||||||||||||||||||||||||||||||||||||||
## BRs, normalized to the SM ones: they scale as (coupling/coupling_SM)^2 / (total/total_SM) self.modelBuilder.factory_('expr::c6_BRscal_hvv("@0*@0/@1", kV, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_htt("@0*@0/@1", ktau, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_hbb("@0*@0/@1", kbottom, c6_Gscal_tot)') self.modelBuilder.factory_('expr::c6_BRscal_hgg("@0*@0/@1", kgamma, c6_Gscal_tot)') if self.doHZg: self.modelBuilder.factory_('expr::c6_BRscal_hzg("@0*@0/@1", kZgamma, c6_Gscal_tot)') | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | def getHiggsSignalYieldScale(self,production,decay,energy): name = "c6_XSBRscal_%s_%s" % (production,decay) print '[LOFullParametrization::C6]' print name, production, decay, energy if self.modelBuilder.out.function(name) = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production = "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay == "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | def getHiggsSignalYieldScale(self,production,decay,energy): name = "c6_XSBRscal_%s_%s" % (production,decay) print '[LOFullParametrization::C6]' print name, production, decay, energy if self.modelBuilder.out.function(name) = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay = "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() | |||||||||||||||||||||||||||||||||||||||||||||
## This lines below are instead in HiggsCouplings.py from LOFullParametrization import C5, C6, C7, PartialWidthsModel c6 = C6()
%ENDCODE%
Working example 1: κγ | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | We will take the full combination but excluding the H→inv search, and extract the κγ modifier, profiling the others. However, to profile them correctly we need to select ranges that are large enough so that the parameter doesn't end on the boundary. We can look at the past CMS Higgs combination, HIG-13-005, to see what was the sensitivity to the other parameters, e.g from the summary plot![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | We will take the full combination but excluding the H→inv search, and extract the κγ modifier, profiling the others. However, to profile them correctly we need to select ranges that are large enough so that the parameter doesn't end on the boundary. We can look at the past CMS Higgs combination, HIG-13-005, to see what was the sensitivity to the other parameters, e.g from the summary plot![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
In addition to the observed outcome, we will compute also the Asimov one | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
text2workspace.py comb.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:c6 -o c6_comb.root --for-fit RANGES="--setPhysicsModelParameterRanges kbottom=0,5:ktop=0,8:ktau=0,4:kgluon=0,4:kgamma=0,3" | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1472 to 1499 | ||||||||||||||||||||||||||||||||||||||||||||||
Working example 2: κg vs κb degeneracy | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | The channels with the best sensitivity to the signal strength are mostly driven by the gluon fusion production cross section (e.g. H→WW and H→ZZ, H→γγboosted H→ττ). Conversely, the H→bb is the dominant decay mode, so that the total width is mostly sensitive to κb. The combination of these two features results in a strong correlation between the uncertainties of the two measurements of κg and κb: a simultaneous increase of both couplings would produce an increase in the gluon fusion production cross section but also a decrease of the branching ratio in the other modes (WW, ZZ, γγ and ττ), leaving the expected σ¡¤BR unchanged, and so provides still a good fit to the data. | |||||||||||||||||||||||||||||||||||||||||||||
> > | The channels with the best sensitivity to the signal strength are mostly driven by the gluon fusion production cross section (e.g. H→WW and H→ZZ, H→γγboosted H→ττ). Conversely, the H→bb is the dominant decay mode, so that the total width is mostly sensitive to κb. The combination of these two features results in a strong correlation between the uncertainties of the two measurements of κg and κb: a simultaneous increase of both couplings would produce an increase in the gluon fusion production cross section but also a decrease of the branching ratio in the other modes (WW, ZZ, γγ and ττ), leaving the expected σ·BR unchanged, and so provides still a good fit to the data. | |||||||||||||||||||||||||||||||||||||||||||||
We can see this feature with a 2D likelihood scan: | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
./parallelScan.py c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgluon -P kbottom --floatOtherPOI=1 $RANGES --points=400 -n C6_comb_SCAN_2D_kgluon_kbottom | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1500 to 1529 | ||||||||||||||||||||||||||||||||||||||||||||||
A macro to make the summary plot is below: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_all.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_all() { TCanvas c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_all.root", "RECREATE"); const int decays = 6; const char *decay[decays] = { "hbb", "htt", "hzz", "hww", "hgg", "comb" }; const int color[decays] = { 67 , 221 , 2 , 4 , 209 , 1 }; const int fill[decays] = { 228 , 223 , 208 , 216 , 211 , 203 }; TList *c68[decays] = { 0,0,0,0,0, 0 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",decay[i])); gFile->Close(); } fOut->cd(); fOut->ls(); for (int i = 0; i < decays; ++i) { c68[i] = (TList) fOut->Get(Form("kvkf_%s_c68",decay[i])); styleMultiGraph(c68[i], color[i], 3, 1); fillMultiGraph(c68[i], fill[i], 1001); } TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_all.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_all() { TCanvas c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_all.root", "RECREATE"); const int decays = 6; const char *decay[decays] = { "hbb", "htt", "hzz", "hww", "hgg", "comb" }; const int color[decays] = { 67 , 221 , 2 , 4 , 209 , 1 }; const int fill[decays] = { 228 , 223 , 208 , 216 , 211 , 203 }; TList *c68[decays] = { 0,0,0,0,0, 0 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",decay[i])); gFile->Close(); } fOut->cd(); fOut->ls(); for (int i = 0; i < decays; ++i) { c68[i] = (TList) fOut->Get(Form("kvkf_%s_c68",decay[i])); styleMultiGraph(c68[i], color[i], 3, 1); fillMultiGraph(c68[i], fill[i], 1001); } TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | for (int i = 0; i < decays; ++i) { c68[i]->Draw("F SAME"); } for (int i = 0; i < decays; ++i) { c68[i]->Draw("L SAME"); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | for (int i = 0; i < decays; ++i) { c68[i]->Draw("F SAME"); } for (int i = 0; i < decays; ++i) { c68[i]->Draw("L SAME"); } | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TGraph *best = fOut->Get("kvkf_comb_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TGraph *best = fOut->Get("kvkf_comb_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_kvkf_all.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_kvkf_all.png"); fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1543 to 1573 | ||||||||||||||||||||||||||||||||||||||||||||||
kgluon : +1.039 kgamma : +1.143 Done in 0.11 min (cpu), 0.11 min (real) | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
combine c6_comb.root -M MultiDimFit --algo=singles -m 125.7 $RANGES | |||||||||||||||||||||||||||||||||||||||||||||
> > | combine c6_comb.root -M MultiDimFit --algo=singles -m 125.7 $RANGES | |||||||||||||||||||||||||||||||||||||||||||||
>>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1617 to 1645 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp" }% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"python"}% class C6I(SMLikeHiggsModel): "as C6 but including direct searches for H->invis" def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False self.doHZg = False def setPhysicsOptions(self,physOptions): for po in physOptions: if po.startswith("higgsMassRange="): self.floatMass = True self.mHRange = po.replace("higgsMassRange=","").split(",") print 'The Higgs mass range:', self.mHRange if len(self.mHRange) = 2: raise RuntimeError, "Higgs mass range definition requires two extrema" elif float(self.mHRange[0]) >= float(self.mHRange[1]): raise RuntimeError, "Extrama for Higgs mass range defined with inverterd order. Second must be larger the first" if po = 'doHZg': self.doHZg = True def doParametersOfInterest(self): self.modelBuilder.doVar("kV[1,0.0,2.0]") self.modelBuilder.doVar("ktau[1,0.0,2.0]") self.modelBuilder.doVar("ktop[1,0.0,2.0]") self.modelBuilder.doVar("kbottom[1,0.0,2.0]") self.modelBuilder.doVar("kgluon[1,0.0,2.0]") self.modelBuilder.doVar("kgamma[1,0.0,2.0]") self.modelBuilder.doVar("BRInv[0,0,1]") pois = 'kV,ktau,ktop,kbottom,kgluon,kgamma,BRInv' if self.doHZg: self.modelBuilder.doVar("kZgamma[1,0.0,30.0]") pois + ",kZgamma" if self.floatMass: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setRange(float(self.mHRange[0]),float(self.mHRange[1])) self.modelBuilder.out.var("MH").setConstant(False) else: self.modelBuilder.doVar("MH[%s,%s]" % (self.mHRange[0],self.mHRange[1])) self.modelBuilder.doSet("POI",pois+',MH') else: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setVal(self.options.mass) self.modelBuilder.out.var("MH").setConstant(True) else: self.modelBuilder.doVar("MH[%g]" % self.options.mass) self.modelBuilder.doSet("POI",pois) self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup() | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"python"}% class C6I(SMLikeHiggsModel): "as C6 but including direct searches for H->invis" def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False self.doHZg = False def setPhysicsOptions(self,physOptions): for po in physOptions: if po.startswith("higgsMassRange="): self.floatMass = True self.mHRange = po.replace("higgsMassRange=","").split(",") print 'The Higgs mass range:', self.mHRange if len(self.mHRange) = 2: raise RuntimeError, "Higgs mass range definition requires two extrema" elif float(self.mHRange[0]) >= float(self.mHRange[1]): raise RuntimeError, "Extrama for Higgs mass range defined with inverterd order. Second must be larger the first" if po = 'doHZg': self.doHZg = True def doParametersOfInterest(self): self.modelBuilder.doVar("kV[1,0.0,2.0]") self.modelBuilder.doVar("ktau[1,0.0,2.0]") self.modelBuilder.doVar("ktop[1,0.0,2.0]") self.modelBuilder.doVar("kbottom[1,0.0,2.0]") self.modelBuilder.doVar("kgluon[1,0.0,2.0]") self.modelBuilder.doVar("kgamma[1,0.0,2.0]") self.modelBuilder.doVar("BRInv[0,0,1]") pois = 'kV,ktau,ktop,kbottom,kgluon,kgamma,BRInv' if self.doHZg: self.modelBuilder.doVar("kZgamma[1,0.0,30.0]") pois + ",kZgamma" if self.floatMass: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setRange(float(self.mHRange[0]),float(self.mHRange[1])) self.modelBuilder.out.var("MH").setConstant(False) else: self.modelBuilder.doVar("MH[%s,%s]" % (self.mHRange[0],self.mHRange[1])) self.modelBuilder.doSet("POI",pois+',MH') else: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setVal(self.options.mass) self.modelBuilder.out.var("MH").setConstant(True) else: self.modelBuilder.doVar("MH[%g]" % self.options.mass) self.modelBuilder.doSet("POI",pois) self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup() | |||||||||||||||||||||||||||||||||||||||||||||
def setup(self): | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1629 to 1657 | ||||||||||||||||||||||||||||||||||||||||||||||
## The invisible BR is by construction BRInv self.modelBuilder.factory_('expr::c6i_BRscal_hinv("@0", BRInv)') | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | def getHiggsSignalYieldScale(self,production,decay,energy): name = "c6i_XSBRscal_%s_%s" % (production,decay) if self.modelBuilder.out.function(name) = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production = "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt", "hinv"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay == "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | def getHiggsSignalYieldScale(self,production,decay,energy): name = "c6i_XSBRscal_%s_%s" % (production,decay) if self.modelBuilder.out.function(name) = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt", "hinv"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay = "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() | |||||||||||||||||||||||||||||||||||||||||||||
c6i = C6I() %ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1647 to 1675 | ||||||||||||||||||||||||||||||||||||||||||||||
hideimgleft="/pub/TWiki/TWikiDocGraphics/toggleclose-small.gif" class="twikiHelp" }% | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
text2workspace.py hinv_vbf_8TeV.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_vbfinv.root --for-fit combine -m 125.7 -M MultiDimFit c6i_vbfinv.root -n C6I_vbfinv_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=400 --algo=grid3x3 $RANGES | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
text2workspace.py comb_hww.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_hww.root --for-fit | |||||||||||||||||||||||||||||||||||||||||||||
> > | text2workspace.py comb_hww.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_hww.root --for-fit | |||||||||||||||||||||||||||||||||||||||||||||
./parallelScan.py -M MultiDimFit c6i_hww.root -n C6I_hww_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=100 --algo=grid3x3 $RANGES hadd -f higgsCombineC6I_hww_SCAN_2D_kV_BRInv.MultiDimFit.mH125.7.root higgsCombineC6I_hww_SCAN_2D_kV_BRInv.*.MultiDimFit.mH125.7.root | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
HWW="hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt" | |||||||||||||||||||||||||||||||||||||||||||||
> > | HWW="hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt" | |||||||||||||||||||||||||||||||||||||||||||||
HIN="hinv_8TeV=hinv_vbf_8TeV.txt" combineCards.py -S $HWW $HIN > comb_test2.txt text2workspace.py comb_test2.txt -m 125.7 -P LOFullParametrization:c6i -o c6i_test2.root --for-fit | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1665 to 1690 | ||||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin-->
| ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1685 to 1708 | ||||||||||||||||||||||||||||||||||||||||||||||
./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=1 --points=400 --algo=grid $RANGES ---hadd | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() The solid-filled areas are the 68% and 95% CL regions for the combination not including the direct search for H→invis, the hatched areas are the 68% and 95% CL regions including it. You can appreciate how without the constraint there is a full degeneracy between the measured coupling and the invisible BR, which is then broken by the direct search. | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() The solid-filled areas are the 68% and 95% CL regions for the combination not including the direct search for H→invis, the hatched areas are the 68% and 95% CL regions including it. You can appreciate how without the constraint there is a full degeneracy between the measured coupling and the invisible BR, which is then broken by the direct search. | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1756 to 1782 | ||||||||||||||||||||||||||||||||||||||||||||||
Results:
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
%TWISTY{ mode="div" showlink="Plotting macro" | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1769 to 1794 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp" }% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% void plot1ds_c6i() { TCanvas c1 = new TCanvas("c1","c1"); gStyle->SetOptStat(0); TFile *fOut = TFile::Open("plot1d_c6i.root", "RECREATE"); const int params = 7; const char *param[params] = { "kV", "kbottom", "ktop", "ktau", "kgluon", "kgamma", "BRInv" }; const float min[params] = { 0, 0, 0, 0, 0, 0, 0, }; const float max[params] = { 2, 5, 8, 4, 4, 3, 1, }; const char *label[params] = { "#kappa_{V}", "#kappa_{b}", "#kappa_{t}", "#kappa_{#tau}", "#kappa_{g}", "#kappa_{#gamma}","BR_{Inv}" }; for (int i = 0; i < params; ++i) { cout << param[i] << endl; gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s.MultiDimFit.mH125.7.root", param[i])); TTree *limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph gO = gROOT->FindObject("Graph")->Clone(); gO->Sort(); gFile->Close(); gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s_ASIMOV.MultiDimFit.mH125.7.root", param[i])); limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph *gA = gROOT->FindObject("Graph")->Clone(); gA->Sort(); gFile->Close(); | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% void plot1ds_c6i() { TCanvas c1 = new TCanvas("c1","c1"); gStyle->SetOptStat(0); TFile *fOut = TFile::Open("plot1d_c6i.root", "RECREATE"); const int params = 7; const char *param[params] = { "kV", "kbottom", "ktop", "ktau", "kgluon", "kgamma", "BRInv" }; const float min[params] = { 0, 0, 0, 0, 0, 0, 0, }; const float max[params] = { 2, 5, 8, 4, 4, 3, 1, }; const char *label[params] = { "#kappa_{V}", "#kappa_{b}", "#kappa_{t}", "#kappa_{#tau}", "#kappa_{g}", "#kappa_{#gamma}","BR_{Inv}" }; for (int i = 0; i < params; ++i) { cout << param[i] << endl; gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s.MultiDimFit.mH125.7.root", param[i])); TTree *limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph gO = gROOT->FindObject("Graph")->Clone(); gO->Sort(); gFile->Close(); gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s_ASIMOV.MultiDimFit.mH125.7.root", param[i])); limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph *gA = gROOT->FindObject("Graph")->Clone(); gA->Sort(); gFile->Close(); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | gO->SetName(Form("c6i_1d_%s_obs", param[i])); gA->SetName(Form("c6i_1d_%s_asimov",param[i])); fOut->WriteTObject(gO); fOut->WriteTObject(gA); | |||||||||||||||||||||||||||||||||||||||||||||
> > | gO->SetName(Form("c6i_1d_%s_obs", param[i])); gA->SetName(Form("c6i_1d_%s_asimov",param[i])); fOut->WriteTObject(gO); fOut->WriteTObject(gA); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | TH1F *frame = new TH1F("frame",Form(";%s;- 2 #Delta ln L",label[i]), 10, min[i], max[i]); frame->Draw(); frame->GetYaxis()->SetRangeUser(0,7); gO->SetLineWidth(3); gO->SetLineColor(1); gA->SetLineWidth(2); gA->SetLineColor(1); gA->SetLineStyle(2); gO->Draw("C SAME"); gA->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
> > | TH1F *frame = new TH1F("frame",Form(";%s;- 2 #Delta ln L",label[i]), 10, min[i], max[i]); frame->Draw(); frame->GetYaxis()->SetRangeUser(0,7); gO->SetLineWidth(3); gO->SetLineColor(1); gA->SetLineWidth(2); gA->SetLineColor(1); gA->SetLineStyle(2); gO->Draw("C SAME"); gA->Draw("C SAME"); | |||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | c1->Print(Form("output_C6I_1D_%s.png", param[i])); delete frame; } fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | c1->Print(Form("output_C6I_1D_%s.png", param[i])); delete frame; } fOut->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1794 to 1819 | ||||||||||||||||||||||||||||||||||||||||||||||
The tdefinition of transverse mass mT(ll,ETmiss) used in this analysis is mT(ll,ETmiss) := sqrt( 2 * pT(ll) * ETmiss * (1 - cos Δ#phi;(ll, ETmiss) )
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
Exercise for the reader | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1880 to 1905 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% //// == rootlogon === { TString makeshared(gSystem->GetMakeSharedLib()); TString dummy = makeshared.ReplaceAll("-W ", ""); dummy.ReplaceAll("g++ ","g++ -std=c++0x "); dummy.ReplaceAll("-Wshadow "," "); dummy.ReplaceAll("-Wall ","-Wall -Wno-deprecated-declarations "); gSystem->SetMakeSharedLib(dummy); gSystem->Load("libGenVector.so"); // strangely this is needed if (getenv("CMSSW_BASE")) { gSystem->Load("libFWCoreFWLite.so"); gSystem->Load("libDataFormatsFWLite.so"); gSystem->SetIncludePath("-I$ROOFITSYS/include"); AutoLibraryLoader::enable(); } } %ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% //// = makePlots_PAT.cxx++ = #include #include #include #include #include #include #include #include #include #include #include | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% //// = makePlots_PAT.cxx++ = #include #include #include #include #include #include #include #include #include #include #include | |||||||||||||||||||||||||||||||||||||||||||||
// hide all from ROOT #if defined( CINT) && defined( MAKECINT) #include | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1914 to 1939 | ||||||||||||||||||||||||||||||||||||||||||||||
stackPlots(hmll_sm, hmll_2p, "mll"); stackPlots(hdphill_sm, hdphill_2p, "dphill"); stackPlots(hmT_sm, hmT_2p, "mT"); stackPlots(hm2d_sm, hm2d_2p, "m2d"); | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | output->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
> > | output->Close(); } | |||||||||||||||||||||||||||||||||||||||||||||
/* Note: you can also compile it with gcc, though it's a bit a pain in the neck gcc -std=c++0x -Wno-deprecated -Wno-deprecated-declarations -lstdc++ $(root-config --cflags --ldflags --libs) -lGenVector -I $CMSSW_RELEASE_BASE/src -L $CMSSW_RELEASE_BASE/lib/$SCRAM_ARCH -lDataFormatsPatCandidates -lDataFormatsFWLite -lDataFormatsCommon -lFWCoreUtilities -lFWCoreFWLite -lDataFormatsMath -I$(scramv1 tool tag clhepheader INCLUDE) makePlots_PAT.cxx ; ./a.out */ int main(int argc, char **argv) { makePlots_PAT(); } %ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1934 to 1959 | ||||||||||||||||||||||||||||||||||||||||||||||
You can find the example of such distributions below for the ZZ+WW combination performed in HIG-13-005![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Measures of separation and the CLs criterium | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1953 to 1978 | ||||||||||||||||||||||||||||||||||||||||||||||
![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
It can be seen that the when going from left to right, the distributions P(q|A) and P(q|N) become closer, less discernible, and the separation is smaller. This is due to the fact that the WW channel has a deficit of signal.
Finally, in the following table you can find the results of the hypothesis test performed for HIG-13-005![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Exercise for the reader | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1974 to 1999 | ||||||||||||||||||||||||||||||||||||||||||||||
class="twikiHelp" }% The answer is deceivingly simple: since P( q < qobs | N ) < 0.5, one needs to have signed deviations to note that qobs is to the left of the median of P( q | N ). In this case, since P( q < qobs | N ) < 0.5 and P( q > qobs | A ) > 0.5, it means that the observation is somewhere in between the qexp(A) and qexp(N) as can be seen in the following figure: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
which is part of the plots for the individual ZZ and WW channels![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 1982 to 2007 | ||||||||||||||||||||||||||||||||||||||||||||||
Practice | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | For the spin-2 (2+m in this case) vs spin-0 (0+) hypotheses test, the full H→WW→2a„¡°2ν analysis includes 4 channels corresponding to {0-jet,1-jet}⊗{7 TeV, 8 TeV}. For simplicity we will perform the exercise on the 0-jet 8 TeV data card. | |||||||||||||||||||||||||||||||||||||||||||||
> > | For the spin-2 (2+m in this case) vs spin-0 (0+) hypotheses test, the full H→WW→2a„“2ν analysis includes 4 channels corresponding to {0-jet,1-jet}⊗{7 TeV, 8 TeV}. For simplicity we will perform the exercise on the 0-jet 8 TeV data card. | |||||||||||||||||||||||||||||||||||||||||||||
head jcp_hwwof_0j_8TeV.txt [...] Observation 5747 | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2004 to 2027 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% //Compare main signal and background | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | histo_ggH->GetXaxis()->SetRangeUser(-1,1) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_qqWW->SetLineColor(kGray) histo_qqWW->DrawNormalized("Hsame") | |||||||||||||||||||||||||||||||||||||||||||||
> > | histo_ggH->GetXaxis()->SetRangeUser(-1,1) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_qqWW->SetLineColor(kGray) histo_qqWW->DrawNormalized("Hsame") | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
You will notice multiple peaks. These correspond to an unrolled 2D histogram encoding the information of mT and mll as explained in https://twiki.cern.ch/twiki/bin/view/CMSPublic/Hig12042TWiki![]() | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Let's now compare the different signal model shapes: | ||||||||||||||||||||||||||||||||||||||||||||||
Deleted: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ||||||||||||||||||||||||||||||||||||||||||||||
root -l jcp_hwwof_0j.input_8TeV.root | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2026 to 2048 | ||||||||||||||||||||||||||||||||||||||||||||||
%CODE{"cpp"}% //Define a region where it is easy to see the differences between the alternate signals | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | histo_ggH->;GetXaxis()->SetRangeUser(-1,0) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_ggH_ALT->SetLineColor(kBlue) histo_ggH_ALT->DrawNormalized("Hsame") histo_qqbarH_ALT->SetLineColor(kRed) histo_qqbarH_ALT->DrawNormalized("Hsame") | |||||||||||||||||||||||||||||||||||||||||||||
> > | histo_ggH->;GetXaxis()->SetRangeUser(-1,0) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_ggH_ALT->SetLineColor(kBlue) histo_ggH_ALT->DrawNormalized("Hsame") histo_qqbarH_ALT->SetLineColor(kRed) histo_qqbarH_ALT->DrawNormalized("Hsame") | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
You can see that the shapes are different, with different hypotheses having more or less signal in the different "peaks" (which correspond to different mT values), as shown in: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Building a nested model | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2168 to 2187 | ||||||||||||||||||||||||||||||||||||||||||||||
combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P x --floatOtherPOI=1 --algo=grid -n JCP_X_SCAN_Asimov --points=200 -t -1 --expectSignal=1 --setPhysicsModelParameters x=0,r=1 | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() From there you can see that the data prefers the scalar hypothesis (x=0), and in particular the events are even more scalar-like than what expected on average for a scalar (the likelihood approaches zero with a non-zero slope). You can also see that for larger x the a-priori expected curve grows faster for increasing x (i.e. it has a smaller that the a-priori sensitivity is better) but that is in part compensated by the fact that the observed starts faster. For a gaussian likelihood, a difference in 2*Δ log L of about 3 would correspond to a p-value of about 10%, which is in the same ballpark as the CLs number; a precise comparison cannot be made, since the two kinds of tests are different and since we're considering a variable at the boundary of the interval in which it is defined (and so it is not really expected to behave as a gaussian) | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() From there you can see that the data prefers the scalar hypothesis (x=0), and in particular the events are even more scalar-like than what expected on average for a scalar (the likelihood approaches zero with a non-zero slope). You can also see that for larger x the a-priori expected curve grows faster for increasing x (i.e. it has a smaller that the a-priori sensitivity is better) but that is in part compensated by the fact that the observed starts faster. For a gaussian likelihood, a difference in 2*Δ log L of about 3 would correspond to a p-value of about 10%, which is in the same ballpark as the CLs number; a precise comparison cannot be made, since the two kinds of tests are different and since we're considering a variable at the boundary of the interval in which it is defined (and so it is not really expected to behave as a gaussian) | |||||||||||||||||||||||||||||||||||||||||||||
In all the results above, fqq is not floated: it is a constant parameter set at zero. This means that everything we have done up to now tests gg→2+m vs 0+. | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2193 to 2213 | ||||||||||||||||||||||||||||||||||||||||||||||
The output is for the two cases is:
| ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < |
| |||||||||||||||||||||||||||||||||||||||||||||
> > |
| |||||||||||||||||||||||||||||||||||||||||||||
This means that: - a 2+m graviton produced in qqbar annihilation would be easier to distinguish from the SM Higgs compared to a 2+m graviton produced in gluon fusion - also for the qqbar case the post-fit sensitivity is worse than the pre-fit one (that's understandable since we had already seen that the post-fit signal strength for the SM is below unity, and since the 2+m are not that different) - also for the qqbar case, the observed exclusion is stronger than the expectation, meaning that the data is more SM-like than graviton-like (it's the same data, and the expected distributions for the graviton in the two production modes lie on the same side of the SM, both preferring smaller mT and larger m(ll) ) | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | Again, one can see that also from looking at the likelihood as function of x in th two cases: ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | Again, one can see that also from looking at the likelihood as function of x in th two cases: ![]() | |||||||||||||||||||||||||||||||||||||||||||||
</> <--/twistyPlugin--> | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2224 to 2242 | ||||||||||||||||||||||||||||||||||||||||||||||
and let's explore: | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | %CODE{"cpp"}% q->Print() | |||||||||||||||||||||||||||||||||||||||||||||
> > | %CODE{"cpp"}% q->Print() | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > | ||||||||||||||||||||||||||||||||||||||||||||||
****************************************************************************** *Tree :q : Test statistics * *Entries : 1993 : Total = 42889 bytes File Size = 8799 * | ||||||||||||||||||||||||||||||||||||||||||||||
Line: 2256 to 2275 | ||||||||||||||||||||||||||||||||||||||||||||||
The type encode which hypothesis the value of q corresponds to (observation has type zero). To see that, run:
%CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | q->Draw("type") gPad->SetLogy() | |||||||||||||||||||||||||||||||||||||||||||||
> > | q->Draw("type") gPad->SetLogy() | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
And finally, you can quickly visualise the results of the hypothesis test by running: %CODE{"cpp"}% | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | gPad->SetLogy(kFalse) q->Draw("-2*q>>hnull(50,-10,10)","type>0","goff") q->Draw("-2*q>>halt(50,-10,10)", "type<0","goff") q->Draw("-2*q>>hobs(50,-10,10)", "type==0","goff") hnull->SetLineColor(kBlue) halt->SetLineColor(kOrange) hobs->SetLineColor(kBlack) hobs->Scale(1e3) hobs->SetFillColor(kBlack) hnull->Draw() halt->Draw("same") hobs->Draw("same") | |||||||||||||||||||||||||||||||||||||||||||||
> > | gPad->SetLogy(kFalse) q->Draw("-2*q>>hnull(50,-10,10)","type>0","goff") q->Draw("-2*q>>halt(50,-10,10)", "type<0","goff") q->Draw("-2*q>>hobs(50,-10,10)", "type==0","goff") hnull->SetLineColor(kBlue) halt->SetLineColor(kOrange) hobs->SetLineColor(kBlack) hobs->Scale(1e3) hobs->SetFillColor(kBlack) hnull->Draw() halt->Draw("same") hobs->Draw("same") | |||||||||||||||||||||||||||||||||||||||||||||
%ENDCODE%
Note how the variable q needs to be multiplied by -2 in order to obtain the test-statistic -2*ln(L(A)/L(N)). | ||||||||||||||||||||||||||||||||||||||||||||||
Changed: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
> > | ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Deleted: | ||||||||||||||||||||||||||||||||||||||||||||||
< < | -- ![]() | |||||||||||||||||||||||||||||||||||||||||||||
Comments | ||||||||||||||||||||||||||||||||||||||||||||||
Added: | ||||||||||||||||||||||||||||||||||||||||||||||
> > |
Line: 1 to 1 | ||||||||
---|---|---|---|---|---|---|---|---|
|
Line: 1 to 1 | ||||||||
---|---|---|---|---|---|---|---|---|
|
Line: 1 to 1 | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Added: | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
> > |
<-- WE DEFINE HERE THE PATH TO THE GITHUB, USED FOR LINKS LATER
Higgs Combination and Properties: Long Exercise - CMSDAS 2014
IntroductionOverviewThe purpose of this long exercise is to understand how some of the key properties of the Higgs boson are measured at CMS: how the couplings are determined from the combination of the results of the searches in different final states (decay modes, production topologies), and how exotic JCP hypotheses are tested and excluded from the kinematics of the Higgs boson decay products. For this exercise, we will use the data from the real higgs searches in CMS, and analyze them with the standard CMS statistical toolkit that has been developed for this purpose. In order to keep the amount of computing time down to a manageable level, we will use only the data from the most sensitive channels from each analysis (but including also new channels that were not present at the time of the latest CMS combination), and a version of the statistical tools containing the newest optimized code that is not yet in the production release.Prerequisites
Course Structure Outline
Miscellaneous NotesThe color scheme of the Exercise is as follows:
Intro I: The Higgs at CMSMeasurements of the Higgs properties in CMS exploit the different production and decay modes expected for the SMH boson. The SMH boson is produced through gluon-fusion (ggH, red), vector-boson fusion (VBF, blue), associated vector-boson production (VH, purple - where the Higgs is radiated off an energetic V=W,Z boson), and in associated production with a top-quark pair (ttH, green):![]() ![]() ![]() ![]() ![]() ![]()
![]() ![]() ![]() Intro 2: DatacardsSetting up your computing environmentThe Higgs statistical codeThe statistical toolkit for the Higgs analyses is provided as a package that is compiled within CMSSW, even if it does not depend on the rest of the CMSSW packages. Tthe code relies critically on the RooFit and RooStats packages, and at the moment it is compatible with the versions packaged in the 6.1.1 and 6.1.2 releases (older releases have bugs, newer releases are backwards-incompatible, not tested and not supported) The code can be retrieved from the github code hosting side. The code for the development version that we'll be using for this CMSDAS school is visible at https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/tree/CMSDAS-2014-CERN![]() scramv1 project CMSSW CMSSW_6_1_2 cd CMSSW_6_1_2/src cmsenv git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit/ cd HiggsAnalysis/CombinedLimit/ git checkout CMSDAS-2014-CERN scramv1 b -j 5after the first installation you can instead just do cd CMSSW_6_1_2/src/HiggsAnalysis/CombinedLimit cmsenv The datacardsThe datacards contain CMS data, so they cannot be put on github. Normally, they're kept in a protected svn repository at CERN, but for convenience you will be able to get them instead from a tarball on AFS.cd CMSSW_6_1_2/src/HiggsAnalysis/CombinedLimit mkdir datacards cd datacards tar xvf /afs/cern.ch/user/g/gpetrucc/public/CMSDAS-2014-CERN-cards.tar.gz A counting experiment: search for invisible higgs decays in the VBF modeA counting experiment is a search where we just count the number of events that pass a selection, and compare that number with the event yields expected from the signal and background. For this kind of analysis, all the information needed for the statistical interpretation of the results is encoded in a simple text file An example of search that is done this way is the one for a Higgs boson produced in the VBF mode and decaying to invisible particles (something that does not happen in the standard model but can happen in other scenarios like supersymmetry) You can look at the datacardhinv_vbf_8TeV.txt
Here is a line-by-line commentary of it
# Invisible Higgs analysis for mH=125 GeVThis line is just a comment, describing the datacard. All lines with "#" are comments, and ignored in the processing. Similarly, lines all of minus characters are ignored imax 1 number of channels jmax 6 number of backgrounds kmax 9 number of nuisance parameters (sources of systematical uncertainties)These lines describe the basics of the datacard: the number of channels, physical processes, and systematical uncertainties Only the first two words (e.g. imax 6 ) are interpreted, the rest is a comment.
------------ # we have just one channel, in which we observe 0 events bin vbf observation 390These lines describe the number of observed events in each final state. In this case, there's a final state, with label vbf that contains 390 events (the comment above "we observe 0 events" is probably a leftover from an old cut-and-paste...)
------------ bin vbf vbf vbf vbf vbf vbf vbf process qqH zvv wmu wel wtau qcd others process 0 1 2 3 4 5 6 rate 207.589 101.879 67.2274 68.2064 54.4137 36.8432 10.3709These lines describe the number of expected events from each physical process in each final state The first two rows of each column identify the final state (always "vbf" in this datacard, since it's the only channel considered) and the physics process. Higgs boson signals in the different production modes have special names that combine can recognise: ggH for gluon fusion, qqH for VBF, WH and ZH for associated production with a weak boson (some older datacards that don't separate the two contain VH for the sum of the two, though this is discouraged), and ttH for associated production with a top anti-top pair. Background processes can have any name.
The third line tells the code if a given process is by default to be interpreted as a signal (process number 0 or negative) or as a background (positive process number). For compatibility with other tools that use the same datacard format, in each final state the process numbers should be different for the different process, and increasing from left to right along the columns.
The fourth line contains the expected event yield. In this case, about 200 signal events would be expected (leftmost column) over a background of about 340 events (sum of all the others).
------------ lumi_8TeV lnN 1.04 - - - - - 1.04 CMS_vbfhinv_qqh_norm lnN 1.12981 - - - - - - CMS_vbfhinv_zvv_stat gmN 12 - 8.4899 - - - - - CMS_vbfhinv_zvv_norm lnN - 1.25331 - - - - - CMS_vbfhinv_wmu_norm lnN - - 1.23748 - - - - CMS_vbfhinv_wel_norm lnN - - - 1.29736 - - - CMS_vbfhinv_wtau_norm lnN - - - - 1.44471 - - CMS_vbfhinv_qcd_norm lnN - - - - - 1.84373 - CMS_vbfhinv_mc_norm lnN - - - - - - 1.30341The rest of the datacard describes the systematical uncertainties, and will be explained more in detail below Systematical uncertaintiesThe way systematical uncertainties are implemented in the Higgs statistics package is by identifying each independent source of uncertainty, and describing the effect it has on the event yields of the different processes. Each source is identified by a name, and in the statistical model it will be associated with a nuisance parameter. An indidual source of uncertainty can have an effect on multiple processes, also across different channels, and all these effects will be correlated (e.g. the uncertainty on the production cross section for Higgs from gluon fusion will affect the expected event yield for this process in all datacards). The size of the effect doesn't have to be the same, e.g. a 1% uncertainty on the muon identification efficiency will have a 2% effect on H→WW→2μ2ν but a 4% effect on H→ZZ→4l, and anti-correlations are also possible (e.g. an increase in b-tagging efficiency will simultaneously and coherently increase the signal yield in final states that require tagged b-jets and decrease the signal yield in final states that require no tagged b-jets). The use of names for each source of uncertainty allow the code to be able to combine multiple datacards recognizing which uncertainties are in common and which are instead separate. The most common model used for systematical uncertainties is the log-normal![]() lnN keyword in combine. The distribution is characterized by a parameter κ, and affects the expected event yields in a multiplicative way: a positive deviation of +1σ correspond to a yield scaled by a factor κ compared to the nominal one, while a negative deviation of -1σ correspond to a scaling by a factor 1/κ. For small uncertainties, the log-normal is approximately a gaussian. If δx/x is the relative uncertainty on the yield, κ can be set to 1+δx/x.
For example, the first two uncertaintes of the VBF H→inv datacard are log-normals
bin vbf vbf vbf vbf vbf vbf vbf process qqH zvv wmu wel wtau qcd others process 0 1 2 3 4 5 6 rate 207.589 101.879 67.2274 68.2064 54.4137 36.8432 10.3709 ------------ lumi_8TeV lnN 1.04 - - - - - 1.04 CMS_vbfhinv_qqh_norm lnN 1.12981 - - - - - -The first one, associated to the normalization of the luminosity in the 8 TeV run ( lumi_8TeV ) has a 4% effect on the expected event yield of the signal and of the "others" background (probably this background is estimated directly from simulations, not normalized to data in control regions, and so the expected yield scales with the luminosity). The dash ( - ) in the columns of the remaining backgrounds signifies that the uncertainty has no effect on those processes.
The second one, CMS_vbfhinv_qqh_norm is an uncertainty associated to the normalization of only the qqH process, of about 13%.
You can note that the name of the second uncertainty has a prefix CMS_vbfhinv_ , while the first one does not: this is used to indicate that the second uncertainty is something pertinent only to the CMS VBF H→inv analysis, while the first uncertainty is common to all CMS or LHC. While the names of the uncertainties don't matter when considering an individual datacard, these naming conventions are very important when combining datacards, to make sure that the effect of an individual source of uncertainty can be tracked across different channels (e.g. other 8 TeV analyses will likely have a lumi_8TeV uncertainty affecting their signal normalization) and that diffent independent sources of uncertainty are not identified in a single one by mistake (which could happen if different datacards used the same name to point to different things).
Another class of systematical uncertainties on the predictions in the signal region are those associated to statistical uncertainties in the control regions used to estimate the backgrounds, or with the limited size of the MC samples used. If the uncertainties are small, Poisson fluctuations become approximately Gaussian, and they be conveniently approximated by log-normal uncertainties with ![]() bin vbf vbf vbf vbf vbf vbf vbf process qqH zvv wmu wel wtau qcd others process 0 1 2 3 4 5 6 rate 207.589 101.879 67.2274 68.2064 54.4137 36.8432 10.3709 ------------ [...] CMS_vbfhinv_zvv_stat gmN 12 - 8.4899 - - - - -The gamma distribution, in the form we use it, is characterized by two parameters: N, the number of events in the control region (or in MC), and α the weight by which each event is multiplied to determine the expected yield in the signal region (if the weights are not all the same, one can take α to be the average weight). By definition, the prediciton in the signal region is gmN immediately followed by the number of events N, and then putting the factor α in the column corresponding to the process. So, in this case N=12, α=8.5 and the expected yield is about 102 events.
Exercise for the reader
What is the total expected background yield, with its statistical and systematical uncertainties?For this approximate computation you can treat the log-normals and Gamma distributions as if they were Gaussians (remembering that approximately κ=1+δx/x and δ x/x ~ 1/sqrt(N+1)). <--/twistyPlugin twikiMakeVisibleInline--> The total yield for the background is 101.879 + 67.2274 + 68.2064 + 54.4137 + 36.8432 + 10.3709 = 338.94 events
The statistical uncertainty is the usual Poisson one, so sqrt(338.94) = 18.4
For the systematical uncertainties, each row is independent, so we can compute the total background uncertainty by adding up the individual contributions in quadrature:
<--/twistyPlugin--> A shape analysis using templates: H→ττA shape analysis relies not only on the expected event yields for the signals and backgrounds, but also on the distribution of those events in some discriminating variable. This approach is often convenient with respect to a counting experiment, especially when the background cannot be predicted reliably a-priori, since the information from the shape allows a better discrimination between a signal-like and a background-like excess, and provides an in-situ normalization of the background. The simpler case of shape analysis is the binned one: the observed distribution of the events in data and the expectations from the signal and all backgrounds are provided as histograms, all with the same binning. Mathematically, this is equivalent to just making N counting experiments, one for each bin of the histogram; however, in practice it's much more convenient to provide the the predictions as histograms directly. For this purpose we will look at the H→ττ search in the τμτh final state with a VBF di-jet topology. The datacard ishtt_mt_5_8TeV.txt , and the full content is:
max 1 number of categories jmax 9 number of samples minus one kmax * number of nuisance parameters ------------------------------------------------------------------------------- shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC shapes ggH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC shapes qqH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC shapes VH * htt_mt.input_8TeV.root $CHANNEL/$PROCESS$MASS $CHANNEL/$PROCESS$MASS_$SYSTEMATIC ------------------------------------------------------------------------------- bin muTau_vbf observation 199 ------------------------------------------------------------------------------- bin muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf process -2 -1 0 1 2 3 4 5 6 7 process ggH qqH VH ZTT QCD W ZL ZJ TT VV rate 2.36529 7.04351 0.028046 94.286 39.2563 46.6949 0.828452 3.72816 4.29862 1.48636 ------------------------------------------------------------------------------- lumi_8TeV lnN 1.042 1.042 1.042 - - - - - - - pdf_qqbar lnN - 1.03 1.02 - - - - - - - pdf_gg lnN 1.08 - - - - - - - - - UEPS lnN 1.04 1.04 1.04 - - - - - - - QCDscale_ggH2in lnN 1.3 - - - - - - - - - QCDscale_qqH lnN - 1.04 - - - - - - - - QCDscale_VH lnN - - 1.04 - - - - - - - CMS_eff_m lnN 1.02 1.02 1.02 1.020 - - 1.020 1.020 1.020 1.020 CMS_eff_t_mutau_8TeV lnN 1.08 1.08 1.08 1.080 - - - - 1.080 1.080 CMS_eff_t_mutau_high_8TeV lnN 1.015 1.015 1.015 1.015 - - - - 1.015 1.015 CMS_scale_t_mutau_8TeV shape 1 1 1 1.000 - - - - - - CMS_scale_j_8TeV lnN 1.1 1.05 1.15 - - - - - 1.050 1.100 CMS_htt_scale_met_8TeV lnN 1.05 1.05 1.05 - - - 1.050 1.050 1.070 1.070 CMS_htt_zttNorm_8TeV lnN - - - 1.030 - - 1.030 1.030 - - CMS_htt_ttbarNorm_8TeV lnN - - - - - - - - 1.080 - CMS_htt_ttbarNorm_mutau_vbf_8TeV lnN - - - - - - - - 1.330 - CMS_htt_DiBosonNorm_8TeV lnN - - - - - - - - - 1.150 CMS_htt_DiBosonNorm_mutau_vbf_8TeV lnN - - - - - - - - - 1.280 CMS_htt_WNorm_mutau_vbf_8TeV lnN - - - - - 1.124 - - - - CMS_htt_extrap_ztt_mutau_vbf_8TeV lnN - - - 1.100 - - - - - - CMS_htt_QCDSyst_mutau_vbf_8TeV lnN - - - - 1.190 - - - - - CMS_htt_ZLeptonFakeTau_mutau_8TeV lnN - - - - - - 1.300 - - - CMS_htt_ZJetFakeTau_mutau_vbf_8TeV lnN - - - - - - - 1.400 - -Compared to the H→inv datacard, you can see that the format is pretty similar, except for two aspects:
shapes lines are used to associate to each channel and physical process the histogram corresponding to the expected distribution of the events, and similarly to associate to the observed data the observed distribution. The format is:
*shapes* _process_ _channel_ _file_ _histogram-name_ _histogram-name-for-systematics_The first two columns after shapes identify to which process and channel the line refers to; the wildcard * can be used to imply that the rule applies to all processes or to all channels, unless a more specific rule exists for it. The channel names are those used in the bin line in the datacard, and the process names are those in the process line, or data_obs for the observed data.
The third and fourth column are used to point to the shape. The patterns $CHANNEL , $PROCESS , $MASS in the file name or histogram name will be substituted by the name of the channel, the name of the process, and the higgs boson mass hypothesis. The fifth column has a similar use, but is used to find the alternative shapes used to encode systematical uncertainties, and will be described later.
Let's now open the rootfile associated with this datacard in root, and look at the histogram associated to the QCD background. There are no specific shapes rules for this QCD background, so the the first, more generic, rule will apply
shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATICand for us $CHANNEL is muTau_vbf and $PROCESS is QCD , so we expect the histogram to be in a directory muTau_vbf and be called QCD .
%CODE{"cpp"}%
/// root -l root -l htt_mt.input_8TeV.root
/// Attaching file htt_mt.input_8TeV.root as _file0...
_file0->cd("muTau_vbf"); QCD->Draw("HIST");
%ENDCODE%
![]() muTau_vbf/data_obs
%CODE{"cpp"}%
cout << data_obs->Integral() << endl;
/// ----> 199 data_obs->Draw();
%ENDCODE%
![]() Form("pattern", args) method of ROOT to create the proper string with the mass written inside.
%CODE{"cpp"}%
qqH110->Draw("HIST");
for (int m = 110; m <= 140; m += 10) { TH1F h = (TH1F) gDirectory->Get(Form("qqH%d",m)); h->SetLineColor(206+4*(m-110)/10); h->SetLineWidth(2); h->Draw("HIST SAME"); }
%ENDCODE%
![]() for mH=110, signal mean = 100.5, rms = 15.1 for mH=115, signal mean = 104.6, rms = 16.8 for mH=120, signal mean = 108.8, rms = 17.5 for mH=125, signal mean = 113.2, rms = 17.6 for mH=130, signal mean = 117.3, rms = 18.0 for mH=135, signal mean = 121.3, rms = 19.1 for mH=140, signal mean = 125.5, rms = 21.1 Shape uncertainties using templatesStill using thehtt_mt_5_8TeV.txt datacard, we can now consider a systematical uncertainty that affects not just the normalization but also the shape of the expected distribution for a process. The uncertainty on the tau energy scale is a good example of this. The relevant lines of the datacard are
shapes * * htt_mt.input_8TeV.root $CHANNEL/$PROCESS $CHANNEL/$PROCESS_$SYSTEMATIC [...] ------------------------------------------------------------------------------- bin muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf muTau_vbf process -2 -1 0 1 2 3 4 5 6 7 process ggH qqH VH ZTT QCD W ZL ZJ TT VV rate 2.36529 7.04351 0.028046 94.286 39.2563 46.6949 0.828452 3.72816 4.29862 1.48636 ------------------------------------------------------------------------------- [...] CMS_scale_t_mutau_8TeV shape 1 1 1 1.000 - - - - - -This tells us that there is a shape uncertainty called CMS_scale_t_mutau_8TeV which affects the signals and the ZTT background (this background is determined embedding simulated taus into Z→μμ events, so any data/mc discrepancy on the energy response for taus will affect the prediction for it just like it does for the simulated taus in signal MC)
For shape analyses using templates, the effect of a shape-affecting systematical uncertainty is encoded in the datacard by providing two additional templates corresponding to ¡ÀNσ shifts of the corresponding parameter. For each process affected by this uncertainty, the column in the datacard will contain a value 1/N. In this case, as often, ¡À1σ shapes have been used, and so the columns of the datacard contain a value of 1.
The associated shapes are determined using the last column of the shapes line, which for this datacard reads $CHANNEL/$PROCESS_$SYSTEMATIC , where $SYSTEMATIC is taken to be the name of the uncertainty in the card, and the postixes Up and Down are added to get the up and down templates. Let's see how much does the mean of the expected ZTT distribution change for the ¡À1σ shifts, and also look at the shapes in ROOT. We can also see that the tau energy scale does not affect only the shape of the background, but also the normalization (probably because if taus are reconstructed as more energetic, they have a higher probability of passing the pT cuts or other kinematic cuts of the analsysis)
%CODE{"cpp"}%
/// root -l htt_mt.input_8TeV.root
/// Attaching file htt_mt.input_8TeV.root as _file0...
_file0->cd("muTau_vbf");
TH1F nominal = (TH1F) gDirectory->Get("ZTT");
TH1F tau_up = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVUp");
TH1F tau_dn = (TH1F) gDirectory->Get("ZTT_CMS_scale_t_mutau_8TeVDown");
printf("Nominal shape: mean = %6.2f (norm: %7.3f events)\n", nominal->GetMean(), nominal->Integral()); // --> Nominal shape: mean = 88.03 (norm: 94.286 events) printf("Up shape: mean = %6.2f (norm: %7.3f events)\n", tau_up->GetMean(), tau_up->Integral()); // --> Up shape: mean = 88.86 (norm: 97.172 events) printf("Down shape: mean = %6.2f (norm: %7.3f events)\n", tau_dn->GetMean(), tau_dn->Integral()); // --> Down shape: mean = 87.20 (norm: 91.995 events)
nominal->SetLineWidth(2); nominal->SetLineColor(kBlack); tau_up->SetLineWidth(2); tau_up->SetLineColor(kBlue); tau_dn->SetLineWidth(2); tau_dn->SetLineColor(kRed); nominal->Draw("HIST"); tau_up->Draw("HIST SAME"); tau_dn->Draw("HIST SAME"); nominal->GetYaxis()->SetRangeUser(0,45);
%ENDCODE%
![]() A parametric shape analysis: H→γγIn some cases, it can be convenient to describe the expected signal and background shapes in terms of analytical functions rather than templates; a typical example are the searches where the signal is apparent as a narrow peak over a smooth continuum background. In this context, uncertainties affecting the shapes of the signal and backgrounds can be implemented naturally as uncertainties on the parameters of those analytical functions. It is also possible to adapt an agnostic approach in which the parameters of the background model are left freely floating in the fit to the data, i.e. only requiring the background to be well described by a smooth function. Technically, this is implemented by means of the RooFit package, that allows writing generic probability density functions, and saving them into ROOT files. The pdfs can be either taken from RooFit's standard library of functions (e.g. Gaussians, polynomials, ...) or hand-coded in C++, and combined together to form even more complex shapes. A prototypical case for this kind of analysis is H→γγ analysis. For this excercise, we will use a datacard that contains only the 8 TeV data for four event categories: cat0 and cat1 are categories of untagged events containing good quality diphotons, the former of the two with a higer purity but lower event yield obtained preferentially selecting high pT diphotons; cat4 and cat5 are categories of di-jet events with different levels of tightness (and so also different level of signal contamination from gluon fusion). The datacard is the following:imax 4 number of bins jmax 5 number of processes minus 1 kmax * number of nuisance parameters ---------------------------------------------------------------------------------------------------------------------------------- shapes WH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_wh_cat0 shapes ZH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_zh_cat0 shapes bkg_mass cat0 hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:pdf_data_pol_model_8TeV_cat0 shapes data_obs cat0 hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:roohist_data_mass_cat0 shapes ggH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_ggh_cat0 shapes qqH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_vbf_cat0 shapes ttH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_tth_cat0 [... same as above for cat1, cat4, cat5 ...] ---------------------------------------------------------------------------------------------------------------------------------- bin cat0 cat1 cat4 cat5 observation -1.0 -1.0 -1.0 -1.0 ---------------------------------------------------------------------------------------------------------------------------------- bin cat0 cat0 cat0 cat0 cat0 cat0 cat1 cat1 cat1 cat1 cat1 cat1 cat4 cat4 cat4 cat4 cat4 cat4 cat5 cat5 cat5 cat5 cat5 cat5 process ZH qqH WH ttH ggH bkg_mass ZH qqH WH ttH ggH bkg_mass ZH qqH WH ttH ggH bkg_mass ZH qqH WH ttH ggH bkg_mass process -4 -3 -2 -1 0 1 -4 -3 -2 -1 0 1 -4 -3 -2 -1 0 1 -4 -3 -2 -1 0 1 rate 6867.0000 19620.0000 12753.0000 19620.0000 19620.0000 1.0000 7259.4000 19620.0000 12360.6000 19620.0000 19620.0000 1.0000 7063.2000 19620.0000 12556.8000 19620.0000 19620.0000 1.0000 3924.0000 19620.0000 15696.0000 19620.0000 19620.0000 1.0000 ---------------------------------------------------------------------------------------------------------------------------------- CMS_eff_j lnN 0.999125 0.964688 0.999125 0.998262 0.996483 - 0.999616 0.980982 0.999616 0.99934 0.999012 - 1.02 1.02 1.02 1.02 1.02 - 1.02 1.02 1.02 1.02 1.02 - CMS_hgg_JECmigration lnN - - - - - - - - - - - - 0.853986 0.995971 0.853986 0.846677 0.927283 - 1.025 1.005 1.025 1.025 1.025 - CMS_hgg_UEPSmigration lnN - - - - - - - - - - - - 0.737174 0.991941 0.737174 0.724019 0.86911 - 1.045 1.01 1.045 1.045 1.045 - CMS_hgg_eff_MET lnN - - - - - - - - - - - - - - - - - - - - - - - - CMS_hgg_eff_e lnN - - - - - - - - - - - - - - - - - - - - - - - - CMS_hgg_eff_m lnN - - - - - - - - - - - - - - - - - - - - - - - - CMS_hgg_eff_trig lnN 1.01 1.01 1.01 1.01 1.01 - 1.01 1.01 1.01 1.01 1.01 - 1.01 1.01 1.01 1.01 1.01 - 1.01 1.01 1.01 1.01 1.01 - CMS_hgg_n_id lnN 1.034/0.958 1.039/0.949 1.034/0.958 1.053/0.915 1.035/0.958 - 1.042/0.936 1.038/0.948 1.042/0.936 1.053/0.909 1.038/0.954 - 1.016/0.972 1.022/0.963 1.016/0.972 1.017/0.967 1.019/0.964 - 1.024/0.959 1.030/0.948 1.024/0.959 1.014/0.975 1.023/0.962 - CMS_hgg_n_pdf_1 lnN - 0.998/0.996 - - 1.002/0.998 - - 0.992/0.999 - - 1.001/1.000 - - 0.999/0.998 - - 1.003/0.999 - - 0.996/1.000 - - 1.002/0.999 - CMS_hgg_n_pdf_10 lnN - 1.028/1.004 - - 0.998/0.998 - - 1.051/0.997 - - 1.004/0.999 - - 1.020/1.000 - - 1.000/0.998 - - 1.010/0.999 - - 0.999/1.000 - [... and more rows like this ...] CMS_hgg_n_sc_gf lnN - - - - 0.842/1.123 - - - - - 0.976/1.031 - - - - - 0.858/1.095 - - - - - 0.880/1.083 - CMS_hgg_n_sc_vbf lnN - 0.993/1.010 - - - - - 0.994/0.994 - - - - - 0.997/1.001 - - - - - 0.999/0.996 - - - - CMS_hgg_n_sigmae lnN 0.950/1.099 0.943/1.114 0.950/1.099 0.956/1.089 0.944/1.112 - 0.954/1.093 0.948/1.104 0.954/1.093 0.971/1.057 0.919/1.165 - 0.996/1.006 0.992/1.016 0.996/1.006 0.995/1.010 0.994/1.012 - 0.991/1.019 0.985/1.031 0.991/1.019 0.995/1.010 0.988/1.026 - CMS_id_eff_eb lnN 1.01999 1.020047 1.01999 1.02002 1.020022 - 1.018535 1.019332 1.018535 1.018642 1.019654 - 1.017762 1.017952 1.017762 1.018824 1.01812 - 1.016723 1.016872 1.016723 1.01884 1.017172 - CMS_id_eff_ee lnN 1.000284 1.000137 1.000284 1.000206 1.0002 - 1.004035 1.001979 1.004035 1.003759 1.001148 - 1.006058 1.005556 1.006058 1.003308 1.005127 - 1.008752 1.008351 1.008752 1.003254 1.007586 - JEC lnN 0.995187 0.938204 0.995187 0.990442 0.980656 - 0.997891 0.966719 0.997891 0.996368 0.994567 - 1.11 1.035 1.11 1.11 1.11 - 1.11 1.035 1.11 1.11 1.11 - QCDscale_VH lnN 0.982/1.021 - 0.982/1.021 - - - 0.982/1.021 - 0.982/1.021 - - - 0.982/1.021 - 0.982/1.021 - - - 0.982/1.021 - 0.982/1.021 - - - QCDscale_ggH lnN - - - - 0.918/1.076 - - - - - 0.918/1.076 - - - - - 0.918/1.076 - - - - - 0.918/1.076 - QCDscale_qqH lnN - 0.992/1.003 - - - - - 0.992/1.003 - - - - - 0.992/1.003 - - - - - 0.992/1.003 - - - - QCDscale_ttH lnN - - - 0.906/1.041 - - - - - 0.906/1.041 - - - - - 0.906/1.041 - - - - - 0.906/1.041 - - UEPS lnN 0.988624 0.858751 0.988624 0.977408 0.954278 - 0.995014 0.923929 0.995014 0.991416 0.987158 - 1.26 1.08 1.26 1.26 1.26 - 1.26 1.08 1.26 1.26 1.26 - lumi_8TeV lnN 1.044 1.044 1.044 1.044 1.044 - 1.044 1.044 1.044 1.044 1.044 - 1.044 1.044 1.044 1.044 1.044 - 1.044 1.044 1.044 1.044 1.044 - pdf_gg lnN - - - 0.920/1.080 0.930/1.076 - - - - 0.920/1.080 0.930/1.076 - - - - 0.920/1.080 0.930/1.076 - - - - 0.920/1.080 0.930/1.076 - pdf_qqbar lnN 0.958/1.042 0.972/1.026 0.958/1.042 - - - 0.958/1.042 0.972/1.026 0.958/1.042 - - - 0.958/1.042 0.972/1.026 0.958/1.042 - - - 0.958/1.042 0.972/1.026 0.958/1.042 - - - vtxEff lnN 0.991/1.025 0.989/1.030 0.991/1.025 0.993/1.020 0.989/1.030 - 0.990/1.026 0.990/1.027 0.990/1.026 0.994/1.016 0.984/1.042 - 1.000/0.999 0.997/1.007 1.000/0.999 1.000/1.003 0.998/1.006 - 0.999/1.004 0.998/1.006 0.999/1.004 1.000/1.000 0.997/1.007 - CMS_hgg_nuissancedeltamcat4 param 0.0 0.001458 CMS_hgg_nuissancedeltafracright_8TeV param 1.0 0.002000 CMS_hgg_nuissancedeltamcat1 param 0.0 0.001470 CMS_hgg_nuissancedeltamcat0 param 0.0 0.001530 CMS_hgg_nuissancedeltasmearcat4 param 0.0 0.001122 CMS_hgg_nuissancedeltasmearcat1 param 0.0 0.001167 CMS_hgg_nuissancedeltasmearcat0 param 0.0 0.001230 CMS_hgg_globalscale param 0.0 0.004717The first difference compared to the template datacard is in the shapes line; let's take for example these two lines
shapes ggH cat0 hgg.inputsig_8TeV_MVA.root wsig_8TeV:hggpdfrel_ggh_cat0 shapes data_obs cat0 hgg.inputbkgdata_8TeV_MVA.root cms_hgg_workspace:roohist_data_mass_cat0In the datacard using templates, the column after the file name would have been the name of the histogram. Here, instead, we found two names, separated by a colon ( : ): the first part identifies the name of the RooWorkspace![]() ![]() ![]() ![]() ![]() Parametric signal normalizationThere is also another feature of the H→γγ datacard that should catch your eye quicky: the event yields are remarkably strange: the the signal yield is19620.0000 events for ggH, qqH and ttH, the background yield is 1.0 , and observed event yield is -1.0 .
Let's start with the simple case: an event yield of -1 just instructs text2workspace and combine to take the yield from the corresponding dataset in the input rootfile, avoiding the need of writing it in the text datacard, but also making the datacard less human-readable. Incidentally, this feature can be used also for datacards that use root histograms like the H→ττ example above.
Now, the signal and background yields: 19k signal events from each production mode with just one background event would be really nice to have, but of course it can't be true. The way the H→γγ works is by relying on an additional feature of text2workspace: in addition to providing generic shapes for the signals and backgrounds, it is also possible to provide generic functions that describe the expected signal and background yields. This additional functions are multiplied by the number in the rate column of the datacard to obtain the final expected event yield. This feature allows e.g. to use the very same datacard to describe all possible different Higgs boson mass hypotheses by just parameterizing properly the expected yield as function of the MH variable.
At present, this feature is not indicated by any special line the datacard: simply, whenever a RooAbsPdf is loaded, text2workspace will search also for a RooAbsReal object with the same name but a _norm postfix, and if present it will use it to scale the event yield.
Armed with this new piece of knowledge, we can now determine what is the expected signal yield for ggH in category 0:
%CODE{"cpp"}%
TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root");
wsig_8TeV->var("MH")->setVal(125.7);
RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0");
RooAbsReal *ggH_norm = wsig_8TeV->function("hggpdfrel_ggh_cat0_norm");
cout << ggH_norm->getVal()*19620.0000 << endl;
// --> 12.3902
%ENDCODE%
This approach can also be used to make a background freely floating, by associating to them as normalization term a floating RooRealVar. However, this can be achieved in a more transparent way by putting instead a normalization uncertainty on that background using a flat pdf: e.g. to leave the background floating between 50% and 200% of its input prediction, a lnU systematic can be used with κ=2 ( lnU has a syntax like lnN , but produces a uniform pdf between 1/κ and κ rather than a log-normal; see SWGuideHiggsAnalysisCombinedLimit![]() Shape uncertainties using parametersThe part of the H→γγ datacard related to the systematics starts with many lines of log-normals that should already be familiar, except possibly for the notation with two numbers separated by a slash (e.g. 0.950/1.099). This notation is used for asymmetrical uncertainties: a log-normal with 0.950/1.099 means that at -1σ the yield is scaled down by a factor 0.95, while at +1σ the yield is scaled up by a factor 1.099. The last part of the datacard contains some lines that use a different syntax, e.g.CMS_hgg_globalscale param 0.0 0.004717These lines encode uncertainties on the parameters of the signal and background pdfs. The example line quoted here informs text2workspace that the parameter CMS_hgg_globalscale is to be assigned a Gaussian uncertainty of ¡À0.004717 around its mean value of zero (0.0).
The effect can be visualized from RooFit, e.g.
%CODE{"cpp"}%
TFile *sig = TFile::Open("hgg.inputsig_8TeV_MVA.root");
wsig_8TeV->var("MH")->setVal(125.7);
RooAbsPdf *ggH = wsig_8TeV->pdf("hggpdfrel_ggh_cat0");
// prepare the canvas
RooPlot *plot = wsig_8TeV->var("CMS_hgg_mass")->frame();
// plot nominal pdf
ggH->plotOn(plot, RooFit::LineColor(kBlack));
// plot minus 3 sigma pdf
wsig_8TeV->var("CMS_hgg_globalscale")->setVal(-3*0.004717);
ggH->plotOn(plot, RooFit::LineColor(kBlue));
// plot plus 3 sigma pdf
wsig_8TeV->var("CMS_hgg_globalscale")->setVal(+3*0.004717);
ggH->plotOn(plot, RooFit::LineColor(kRed));
plot->Draw();
%ENDCODE%
![]() Exercises for the reader
<--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"cpp"}%
TFile *fDat = TFile::Open("hgg.inputbkgdata_8TeV_MVA.root");
cout << cms_hgg_workspace->data("roohist_data_mass_cat0")->sumEntries() << endl;
// .... and so on
%ENDCODE%
The result is:
<--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"cpp"}%
TFile *fSig = TFile::Open("hgg.inputsig_8TeV_MVA.root");
wsig_8TeV->var("MH")->setVal(125.7);
cout << wsig_8TeV->function("hggpdfrel_ggh_cat0_norm")->getVal()*19620.0000 << endl;
cout << wsig_8TeV->function("hggpdfrel_vbf_cat0_norm")->getVal()*19620.0000 + wsig_8TeV->function("hggpdfrel_wh_cat0_norm")->getVal()*12753.0000 + wsig_8TeV->function("hggpdfrel_zh_cat0_norm")->getVal()*6867.0000 << endl; // .... and so on
%ENDCODE%
<--/twistyPlugin--> Unbinned, multi-dimensional analysis: H→ZZThe datacard is identical in form to that of H→γγ, i.e. it provides shapes as RooAbsPdf from a RooWorkspace, and implements shape uncertainties through parameters (for some a fixed range is given to the parameter, e.g.CMS_zz4l_bkgMELA is constrained to be in [-3,3] ; this can be useful if for too large deviations of the parameter the shapes start to misbehave)
imax 1 jmax 7 kmax * ------------ shapes * * hzz4l_4muS_8TeV_0.input.root w:$PROCESS ------------ bin a1_0 observation 13 ------------ ## mass window [105.7,140.7] bin a1_0 a1_0 a1_0 a1_0 a1_0 a1_0 a1_0 a1_0 process ggH qqH WH ZH ttH bkg2d_qqzz bkg2d_ggzz bkg2d_zjets process -4 -3 -2 -1 0 1 2 3 rate 1.0000 1.0000 1.0000 1.0000 1.0000 7.2035 0.1360 0.7119 ------------ lumi_8TeV lnN 1.044 1.044 1.044 1.044 1.044 1.044 1.044 - pdf_gg lnN 1.0720 - - - 1.0780 - 1.0710 - pdf_qqbar lnN - 1.0270 1.0350 1.0349 - 1.0342 - - pdf_hzz4l_accept lnN 1.02 1.02 1.02 1.02 1.02 - - - QCDscale_ggH lnN 1.0750 - - - - - - - QCDscale_qqH lnN - 1.0020 - - - - - - QCDscale_VH lnN - - 1.0040 1.0155 - - - - QCDscale_ttH lnN - - - - 1.0655 - - - QCDscale_ggVV lnN - - - - - - 1.2435 - QCDscale_VV lnN - - - - - 1.0285 - - BRhiggs_hzz4l lnN 1.02 1.02 1.02 1.02 1.02 - - - CMS_eff_m lnN 1.043 1.043 1.043 1.043 1.043 1.043 1.043 - CMS_hzz4mu_Zjets lnN - - - - - - - 0.6/1.4 CMS_zz4l_bkgMELA param 0 1 [-3,3] QCDscale_ggH2in lnN - - - - - - - - QCDscale_qqH lnN - - - - - - - - CMS_zz4l_ggH_PToverM_sys param 0 1 [-3,3] CMS_zz4l_qqH_PToverM_sys param 0 1 [-3,3] CMS_zz4l_ttH_PToverM_sys param 0 1 [-3,3] CMS_zz4l_WH_PToverM_sys param 0 1 [-3,3] CMS_zz4l_ZH_PToverM_sys param 0 1 [-3,3] CMS_zz4l_qqZZ_PToverM_sys param 0 1 [-3,3] CMS_zz4l_ggZZ_PToverM_sys param 0 1 [-3,3] CMS_zz4l_ZX_PToverM_sys param 0 1 [-3,3] CMS_zz4l_mean_m_sig param 0.0 1.0 ## CMS_zz4l_mean_m_sig = 0.001 CMS_zz4l_sigma_m_sig param 0.0 0.2 CMS_zz4l_n_sig_1_8 param 0.0 0.01 interf_ggH param 0 1 [-1,1]The interesting difference however is the observed dataset: %CODE{"cpp"}% gSystem->Load("libHiggsAnalysisCombinedLimit.so") gFile = TFile::Open("hzz4l_4muS_8TeV_0.input.root"); RooAbsData *data = w->data("data_obs"); data->Print("") // --> RooDataSet::data_obs[CMS_zz4l_mass,melaLD,CMS_zz4l_PToverM] = 13 entries %ENDCODE% In the Hγγ the dataset was a RooDataHist with one entry per bin, each entry having a weight equal to the number of events in that bin. Here instead it's a RooDataSet, with one entry per event. Let's explore its contents, both building some binned 1D projections and by making a 2D scatter plot. %CODE{"cpp"}% w->var("CMS_zz4l_mass")->Print(""); // --> RooRealVar::CMS_zz4l_mass = 114.157 L(104 - 142) B(200) RooPlot *mass = w->var("CMS_zz4l_mass")->frame(RooFit::Bins(19)); data->plotOn(mass); mass->Draw(); RooPlot *mela = w->var("melaLD")->frame(RooFit::Bins(10)); data->plotOn(mass); mass->Draw(); TGraph points(data->numEntries()); for (int i = 0; i < data->numEntries(); ++i) { const RooArgSet *ev = data->get(i); points->SetPoint(i, ev->getRealValue("CMS_zz4l_mass"), ev->getRealValue("melaLD")); } points->SetMarkerStyle(20); points->GetXaxis()->SetTitle("m_{4l} (GeV)") points->GetYaxis()->SetTitle("K_{D}") .x tdrstyle.cc points->Draw("AP"); %ENDCODE%
Combining datacardsIn order to analyze the outcomes of multiple searches together, it is necessary to combine their datacards into a single one. This is done by thecombineCards.py script, whose syntax is simple:
combineCards.py Name1=card1.txt Name2=card2.txt .... > card.txtIf the input datacards had just one bin each, then the output channels will be called Name1, Name2, and so on. Otherwise, a prefix Name1_ ... Name2_ will be added to the bin labels in each datacard. The supplied bin names Name1, Name2, etc. must themselves conform to valid C++/python identifier syntax. The package containing the datacards also has a simple script mkcomb.sh which will construct different combinations of input datacards.
For example, for the H→WW search, we will combine the 0-jets opposite-flavour channel and the di-jet opposite-flavour channels, which have respectively the best sensitivity to gluon fusion and vbf higgs production. The combineCards.py $HWW > comb_hww.txt command will be expanded in
combineCards.py hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt > comb_hww.txtThe combined datacard looks like the following Combination of hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt imax 2 number of bins jmax 23 number of processes minus 1 kmax 89 number of nuisance parameters ---------------------------------------------------------------------------------------------------------------------------------- shapes * hww0j_8TeV hwwof_0j.input_8TeV.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC shapes data_obs hww0j_8TeV hwwof_0j.input_8TeV.root histo_Data shapes * hww2j_8TeV hww-19.47fb.mH125.of_2j_shape_mll.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC shapes data_obs hww2j_8TeV hww-19.47fb.mH125.of_2j_shape_mll.root histo_Data ---------------------------------------------------------------------------------------------------------------------------------- bin hww0j_8TeV hww2j_8TeV observation 5729.0 24.0 ---------------------------------------------------------------------------------------------------------------------------------- bin hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww0j_8TeV hww2j_8TeV hww2j_8TeV hww2j_8TeV hww2j_8TeV hww2j_8TeV hww2j_8TeV hww2j_8TeV h process ZH qqH WH ggH qqWW WjetsE Wgamma Top WH_SM WjetsM ggWW ZH_SM VV Wg3l qqH_SM Ztt Zjets ggH_SM qqH ggH WW WWewk Vg WJet Top g process -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 -2 0 15 16 17 18 4 7 rate 1.7481 3.0668 5.8955 237.1930 3969.6300 282.7930 115.5780 498.7010 0.0000 331.7970 210.6270 0.0000 132.6090 167.8020 0.0000 45.9960 0.0000 0.0000 5.0290 1.4400 4.0230 1.9080 0.0389 2.4540 14.0300 0 ---------------------------------------------------------------------------------------------------------------------------------- [...] CMS_8TeV_hww_Top_2j_of_2j_extr lnN - - - - - - - - - - - - - - - - - - - - - - - - 1.188 - [...] CMS_hwwof_0j_MVAWHStatBounding_8TeV shape - - 1.0 - - - - - - - - - - - - - - - - - - - - - - - [....] lumi_8TeV lnN 1.044 1.044 1.044 1.044 1.0 - 1.044 - 1.044 - 1.0 1.044 1.044 1.044 1.044 1.044 - 1.044 1.044 1.044 1.044 1.044 1.044 - - 1You can see that the bin and observation rows now contain two columns (one for each channel), and the common lumi_8TeV uncertainty row affects the columns of of processes from the different cards, while e.g. CMS_8TeV_hww_Top_2j_of_2j_extr or CMS_hwwof_0j_MVAWHStatBounding_8TeV have no effect on the columns of the other channel.
In the last command you see that an option -S is passed to the combineCards.py program. This option is needed when combining together counting experiments (such as VBF H→inv) and shape analyses, and tells the code to reintepret the counting experiment as if it were a shape analysis.
Higgs CouplingsThe basicsFrom Datacards to LikelihoodsThe datacards from each search contain the information about the observed events and the expected yields for the various physical processes, backgrounds and signals. For a SM Higgs of fixed mass, there is no free parameter in the theory, and so nothing to fit. To search for deviations from the SM Higgs predictions, we define a model that can predict different expected yields for the signal, and fit the parameters of this model to the data. We will start with the simplest possible model, with only one parameter μ: the predictions for the signal in a given final state in this model will be μ times the yield expected for a SM Higgs (this parameter can be interpreted as a scale factor for all the production cross sections of the signal, so we sometimes write μ = σ/σSM) In the case of a simple counting experiment with no systematical uncertainties, we can write the likelihood function of the observed data as function of the parameter of interest in terms of the observed events ( N), expected signal events for a SM Higgs ( S), and expected background events ( B)Likelihoods in combine and RooStatsThe first step of the statistical analysis is to convert the text datacard into a likelihood function, in the format used by RooStats. For some simple cases, combine does this automatically when given a text file, but it is possible to do it also manually beforehand - and it's necessary to do it manually whenever we want to use a different model than the one with just μ. The conversion from datacard to likelihood is done with thetext2workspace.py program, also contained in the HiggsAnalysis/CombinedLimit package. The syntax is
text2workspace.py datacard.txt -m mass [ -P physicsModel [ --PO physicsModelOptions .... ] ] [other options] -o output.rootThe physicsModel is the reference to the python class that describes what are the parameters of the model (e.g. μ) and how those parameters affect the yields of the signal and background (e.g. that the signal is multiplied by μ while the background not) The most basic Higgs model, defining a single parameter of interest μ is below. As python and roofit don't like greek letters, μ is called r in the code. The code defines two functions, one that declares the parameters, and one that defines how the expected Higgs boson signal yields for each production mode, decay mode and energy depend on those parameters.
StrictSMLikeHiggsModel is defined in the HiggsAnalysis.CombinedLimit.PhysicsModel![]() strictSMLikeHiggs (see the lines at the bottom of the file), so the argument to pass to the -P option will be HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs
Let's use this model to convert the H→WW + 2 jets datacard into a likelihood
text2workspace.py hwwof_2j_shape_8TeV.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs -o hww2j.rootthe beginning of the datacard was ---------------------------------------------------------------------------------------------------- bin of_2j observation 24 shapes * * hww-19.47fb.mH125.of_2j_shape_mll.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC shapes data_obs * hww-19.47fb.mH125.of_2j_shape_mll.root histo_Data ---------------------------------------------------------------------------------------------------- bin of_2j of_2j of_2j of_2j of_2j of_2j of_2j of_2j of_2j of_2j of_2j process ggH qqH ggWW VgS Vg WJet Top WW WWewk VV DYTT process 0 -1 1 2 3 4 5 6 7 8 9 rate 1.44 5.029 0.4295 0.0195 0.0389 2.454 14.03 4.023 1.908 0.4833 1.879 ----------------------------------------------------------------------------------------------------Now we can open the hww2j.root file in ROOT and see its contents. Start a root prompt with the command root.exe -b -l -n and
%CODE{"cpp"}% // Load the libraries needed to interpret the likelihood gSystem->Load("libHiggsAnalysisCombinedLimit.so"); // Open the file gFile = TFile::Open("hww2j.root"); // Retrieve the workspace from the file RooWorkspace *w = (RooWorkspace *) gFile->Get("w"); // Retrieve the ModelConfig from the workspace RooStats::ModelConfig *mc = (RooStats::ModelConfig *) w->genobj("ModelConfig");
// Now we get the list of parameters of the workspace mc->GetParametersOfInterest()->Print("V"); // the output will be // 1) RooRealVar:: r = 1
// Now we get the expected signal yield for ggH and DYTT cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // the output will be 1.44, and 1.879, as expected
// Now we change the value of r to be 2.0... we'd expect the ggH yield to grow to 2.88 and the DYTT to remain the same w->var("r")->setVal(2.0); cout << w->function("n_exp_final_binof_2j_proc_ggH")->getVal() << endl; // ----> 2.88 cout << w->function("n_exp_final_binof_2j_proc_DYTT")->getVal() << endl; // ----> 1.879
// We can also get the number of observed events, that should match the "observation" line in the datacard cout << w->data("data_obs")->sumEntries() << endl; // ----> 24
%ENDCODE%
One-dimensional fits and profile likelihoodNow that we have written the likelihood as function of a parameter μ we can proceed further to see how well different values of μ describe the observed data. The first step is to find the value of μ that describes the data best, i.e. that maximises the likelihood. This value of μ will be denoted astext2workspace.py comb_hww.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs -o comb_hww.root combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_FITThe output of the program will be >>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 Computing limit starting from observation --- MultiDimFit --- best fit parameter values: r : +0.599 Done in 0.01 min (cpu), 0.02 min (real)The output file from combine will contain a tree with a single entry, corresponding to this value. %CODE{"cpp"}% // root.exe -b -l higgsCombineHWW_FIT.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_FIT.MultiDimFit.mH125.7.root as _file0... root [1] limit->Scan("r"); //************************ //* Row * r * //************************ //* 0 * 0.5986483 * //************************ %ENDCODE% Now, to compare different values of μ, we will use the negative log of the profile likelihood function λ
HWW_SCAN_MU
combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_SCAN_MU --algo=grid --points=100 --setPhysicsModelParameterRanges r=0,3This command will be slower than the previous one, since the code now has to maximise the likelihood 101 times: once first to get the best fit μ then 100 times one for each point along the grid. Some of those fits will not converge at the first attempt, so the putput could contain some lines like Minimization did NOT converge, Edm is above max or Covar was made pos def but this is normally ok.
We can inspect the output from ROOT: the output tree in the file higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root will contain one entry for each of the 101 fits,
%CODE{"cpp"}% // root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root
limit->GetEntries() // (const Long64_t)101
limit->Scan("r:2*deltaNLL") //************************************ //* Row * r * 2*deltaNL * //************************************ //* 0 * 0.5986473 * 0 * //* 1 * 0.0149999 * 5.9730005 * //* 2 * 0.0450000 * 5.3444776 * //* 3 * 0.0750000 * 4.7531828 * //* 4 * 0.1049999 * 4.1995768 * //* 5 * 0.1350000 * 3.6845009 *
limit->Draw("2*deltaNLL:r");
%ENDCODE%
![]() ![]() combine -M MaxLikelihoodFit comb_hww.root -m 125.7 -n HWW_FIT_MU --justFit --robustFit=1The two options on the command line instruct combine to do just a fit for μ ( --justFit ), using a slow but robust algorithm ( --robustFit=1 ; without that, for a datacard with many nuisance parameters the calculation of the uncertainties will often fail with error messages like Error in : Minuit2Minimizer::MINOS failed due to invalid function minimum or Error in : Minos error calculation failed for all parameters ). In the output, you can see various values of r being tried, until it finds the right ones
>>> including systematics >>> method used to compute upper limit is MaxLikelihoodFit >>> random number generator seed is 123456 Computing limit starting from observation Searching for crossing at nll = -18788 in the interval 0.598648, 20 r lvl-here lvl-there stepping 2.538783 -11.83519 +0.49447 1.940135 1.180689 -1.63413 -11.83519 0.582041 0.773260 +0.23757 -1.63413 0.174612 0.947873 -0.38172 +0.23757 0.174612 0.825644 +0.07707 -0.38172 0.052384 0.878028 -0.13010 +0.07707 0.052384 0.841359 +0.02136 -0.13010 0.015715 0.857074 -0.04154 +0.02136 0.015715 0.846074 +0.00325 -0.04154 0.004715 0.850788 -0.01552 +0.00325 0.004715 0.847488 -0.00232 -0.01552 0.001414 0.846498 +0.00158 -0.00232 0.000424 0.846922 -0.00008 +0.00158 0.000424 Searching for crossing at nll = -18788 in the interval 0.598648, 0 r lvl-here lvl-there stepping 0.538783 +0.45835 +0.49447 -0.059865 0.478919 +0.35676 +0.45835 -0.059865 0.419054 +0.19666 +0.35676 -0.059865 0.359189 -0.04147 +0.19666 -0.059865 0.401094 +0.13250 -0.04147 -0.017959 0.383135 +0.06353 +0.13250 -0.017959 0.365175 -0.01304 +0.06353 -0.017959 0.377747 +0.04198 -0.01304 -0.005388 0.372359 +0.01921 +0.04198 -0.005388 0.366971 -0.00477 +0.01921 -0.005388 0.370743 +0.01215 -0.00477 -0.001616 0.369127 +0.00497 +0.01215 -0.001616 0.367510 -0.00232 +0.00497 -0.001616 0.368642 +0.00280 -0.00232 -0.000485 0.368157 +0.00061 +0.00280 -0.000485 0.367672 -0.00158 +0.00061 -0.000485 0.368011 -0.00005 -0.00158 -0.000145 --- MaxLikelihoodFit --- Best fit r: 0.598648 -0.230684/+0.248194 (68% CL) nll S+B -> -8.65939 nll B -> -1 Done in 0.11 min (cpu), 0.12 min (real)So, from this combination we have measured Exercise for the reader
Do a fit of the signal strength separately for the 0-jet and 2-jet datacards of H→WW, and compute the weighted average using the textbook formulas for averaging numbers accounting for their uncertainties (for this approximate computation, you can deal with the asymmetric uncertainties by just taking the one that goes in the direction of the average result). How different it is from what obtained from the combined datacard?
<--/twistyPlugin twikiMakeVisibleInline-->
First we convert the datacards in likelihoods
text2workspace.py -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs hwwof_0j_shape_8TeV.txt -o hww0j.root text2workspace.py -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs hwwof_2j_shape_8TeV.txt -o hww2j.rootThen we run the fits on them
(0.636522/(0.267546)^2 + 0.449522/(0.554423)^2)/(1/(0.267546)^2 + 1/(0.554423)^2) = 0.60 with an uncertainty of 1/sqrt(1/(0.267546)^2 + 1/(0.554423)^2) = 0.24, very similar to the result obtained from the fit to the combined datacard. This is understandable, since the systematical uncertainties in common between the two analyses are not the ones driving the sensitivity. <--/twistyPlugin--> Expectations: the Asimov datasetDetermining what is the expected outcome of a fit if the data were to match the SM Higgs expectations can be useful for various purposes, e.g. which to quantify the sensitivity of a given analysis in a way that is independent of the presence of fluctuations in the observed data (or when the data is still blinded, or not yet collected). The behaviour ofa fit with a complex model on the Asimov is also easier to understand with respect to the fit to the observed data, especially if the data departs significantly from the expectations. Normally, the procedure used to determine these expectations is to generate a large sample of possible experimental outcomes, all under the SM Higgs hypothesis, and use those to determine the median outcome of the analysis. This however would require substantial computational resources. Fortunately, a much cheaper but good approximation exists, through what is called the "Asimov dataset" In a parallel reality imagined by the science fiction writer I. Asimov, politics was run in a peculiar way: instead of mobilizing millions of people to cast their vote to deliberate on something, an algorithm was used to select an individual "average" person, and then this person was asked to take the decision on that matter. In more practical terms, an Asimov dataset can be defined as the one that matches perfectly the median expectations from the SM Higgs (i.e. the event yields in each bin are exactly as expected). The idea is simple and intuitive, and apart from a few subtleties like the need of having a non-integer number of observed events in each bin and of dealing efficiently with multi-dimensional models, which we won't discuss in this exercise. When running on a datacard, combine can be instructed to use an Asimov dataset instead of the observed one simply by providing the options-t -1 --expectSignal=1 (the last of which instructs the code to create an Asimov for signal+background and not for background-only).
We can for example produce the expected likelihood scan vs μ for H→WW with
combine -M MultiDimFit comb_hww.root -m 125.7 -n HWW_SCAN_MU_ASIMOV --algo=grid --points=100 --setPhysicsModelParameterRanges r=0,3 -t -1 --expectSignal=1A comparison of the two likelihood scans, observed and Asimov, is below: the observed data is shown as a thick line, the Asimov as a light dotted line. The shapes are similar, but the Asimov is centered at μ=1 (by construction, since that's the SM expectation), while the observed data is not. ![]() <--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"cpp"}% /// root -l higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root // Attaching file higgsCombineHWW_SCAN_MU.MultiDimFit.mH125.7.root as _file0... // Attaching file higgsCombineHWW_SCAN_MU_ASIMOV.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graph = gROOT->FindObject("Graph")->Clone(); graph->Sort();
_file1->cd(); limit->Draw("2*deltaNLL:r"); TGraph *graphA = gROOT->FindObject("Graph")->Clone(); graphA->Sort();
graph->SetLineWidth(3); graph->SetLineColor(1); graphA->SetLineWidth(1); graphA->SetLineColor(4); graphA->SetLineStyle(2); graph->Draw("AC"); graphA->Draw("C SAME");
graph->SetTitle("Likelihood scan as function of #mu;#mu;-2 #Delta ln L"); graph->GetXaxis()->SetRangeUser(0,2); graph->GetYaxis()->SetRangeUser(0,8);
%ENDCODE%
<--/twistyPlugin--> Two-dimensional fits: μggH,ttH vs μVBF,VHNow we will look at a more elaborate model, which has two parameter:
![]() RF in the code) and μVBF,VH ( RV in the code), search for the maximum, and define a text2workspace.py comb_hgg.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:rVrFXSHiggs -o rVrF_hgg.root combine -M MultiDimFit rVrF_hgg.root -m 125.7 -n HGG_RVRF_FIT >>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods. Computing limit starting from observation --- MultiDimFit --- best fit parameter values: RV : +0.259 RF : +1.229 Done in 0.01 min (cpu), 0.01 min (real)Now, for the 2D scan we'll do 400 points (20x20), using the ranges [-1,3] and [-2,5] for the two parameters. This requires 400+1 fits, so it will take a few minutes. combine -M MultiDimFit rVrF_hgg.root -m 125.7 -n HGG_RVRF_SCAN_2D --algo=grid --points=400 --setPhysicsModelParameterRanges RF=-1,3:RV=-2,5 >>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods. Set Range of Parameter RF To : (-1,3) Set Range of Parameter RV To : (-2,5) Computing limit starting from observation Point 0/400, (i,j) = (0,0), RV = -1.825000, RF = -0.900000 Point 2/400, (i,j) = (0,2), RV = -1.825000, RF = -0.500000 Point 4/400, (i,j) = (0,4), RV = -1.825000, RF = -0.100000 Point 6/400, (i,j) = (0,6), RV = -1.825000, RF = 0.300000 Point 8/400, (i,j) = (0,8), RV = -1.825000, RF = 0.700000 Point 10/400, (i,j) = (0,10), RV = -1.825000, RF = 1.100000 Point 12/400, (i,j) = (0,12), RV = -1.825000, RF = 1.500000 [...] Point 388/400, (i,j) = (19,8), RV = 4.825000, RF = 0.700000 Point 390/400, (i,j) = (19,10), RV = 4.825000, RF = 1.100000 Point 392/400, (i,j) = (19,12), RV = 4.825000, RF = 1.500000 Point 394/400, (i,j) = (19,14), RV = 4.825000, RF = 1.900000 Point 396/400, (i,j) = (19,16), RV = 4.825000, RF = 2.300000 Point 398/400, (i,j) = (19,18), RV = 4.825000, RF = 2.700000 Done in 1.52 min (cpu), 1.53 min (real)A simple plot from the output of the bidimensional scan can be done with %CODE{"cpp"}% // plot the likelihood in the 20x20 grid, selecting the interesting region of small deltaNLL limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ"); // make the plot less horrible gStyle->SetOptStat(0); h2->SetContour(200); h2->SetTitle("Profile Likelihood in 2D;#mu_{ggH,ttH};#mu_{VBF,VH}"); /// Now we add the best fit point limit->Draw("RV:RF","deltaNLL == 0","P SAME"); TGraph *bestFit = gROOT->FindObject("Graph"); bestFit->SetMarkerStyle(34); bestFit->SetMarkerSize(3); c1->Update(); %ENDCODE% ![]() ![]() ROOT::Math::chisquared_quantile(cl,ndf) can be used to compute the thresholds for each CL.
ROOT provides a way to extract a smooth contour out of a two-dimensional histogram, and this functionality is packaged in a script test/plotting/contours2D.cxx![]()
contours2D.cxx script in your working directory, and then
%CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root contours2D.cxx // Attaching file higgsCombineHGG_RVRF_SCAN_2D.MultiDimFit.mH125.7.root as _file0... // Processing contours2D.cxx...
// before we had done // limit->Draw("2*deltaNLL:RV:RF>>h2(20,-1,3,20,-2,5)","deltaNLL > 0 && deltaNLL < 11","PROF COLZ");
// now we do
contour2D("RF",20,-1.,3., "RV",20,-2.,5., 1.,1.);
%ENDCODE%
![]() Fitting only a subset of the parameters of a modelSometimes, one is interested in determining confidence regions for one or two parameters in a model that has more parameters. There are two conceptually different ways in which this can be done: the other parameters can be either kept fixed at some nominal value (e.g. the SM prediction), or they can be profiled (i.e. left freely floating in the fit) Still using the H→γγ analysis with the μggH,ttH vs μVBF,VH model as an example, we can compare the two approaches by making 1d likelihood scans as function of μggH,ttH keeping μVBF,VH fixed or profiled: we use the option-P to select the parameter to fit, and the option --floatOtherPOI to configure what to do with the others (μVBF,VH, in our case)
combine -M MultiDimFit rVrF_hgg.root -m 125.7 --algo=grid --points=50 --setPhysicsModelParameterRanges RF=-1,3 -P RF --floatOtherPOI=0 -n HGG_RVRF_SCAN_1D_RF_FIX combine -M MultiDimFit rVrF_hgg.root -m 125.7 --algo=grid --points=50 --setPhysicsModelParameterRanges RF=-1,3 -P RF --floatOtherPOI=1 -n HGG_RVRF_SCAN_1D_RF_PROF%CODE{"cpp"}% // root -l higgsCombineHGG_RVRF_SCAN_1D_RF_* /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_FIX.MultiDimFit.mH125.7.root as _file0... /// Attaching file higgsCombineHGG_RVRF_SCAN_1D_RF_PROF.MultiDimFit.mH125.7.root as _file1... _file0->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphFix = gROOT->FindObject("Graph")->Clone(); graphFix->Sort(); _file1->cd(); limit->Draw("2*deltaNLL:RF"); TGraph *graphProf = gROOT->FindObject("Graph")->Clone(); graphProf->Sort(); graphProf->SetLineWidth(3); graphProf->Draw("AC"); graphFix->SetLineWidth(2); graphFix->SetLineStyle(9); graphFix->SetLineColor(206); graphFix->Draw("C SAME"); %ENDCODE% ![]()
First fits to Higgs search results
<--/twistyPlugin twikiMakeVisibleInline--> First, run the mkcomb.sh to create the combined datacards
bash mkcomb.sh. If you want to get also the combination of the ttH searches, you have to add a line combineCards.py $TTH > comb_ttH.txtThen, for each datacard you can do text2workspace.py -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:strictSMLikeHiggs datacard.txt -o datacard.root combine -M MaxLikelihoodFit --justFit --robustFit=1 -m 125.7 datacard.rootor directly combine -M MaxLikelihoodFit --justFit --robustFit=1 -m 125.7 datacard.txt
<--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline--> for X in h{bb,tt,gg,ww,zz}; do text2workspace.py comb_${X}.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.PhysicsModel:rVrFXSHiggs -o rVrF_${X}.root --for-fit ; done for X in h{bb,tt,gg,ww,zz}; do combine -M MultiDimFit rVrF_${X}.root -m 125.7 -n ${X}_RVRF_SCAN_2D --algo=grid3x3 --points=400 --setPhysicsModelParameterRanges RF=-1,3:RV=-2,5; doneTo see a single plot, interactively in ROOT, %CODE{"cpp"}% /// root -l tdrstyle.cc contours2D.cxx higgsCombinehww_RVRF_SCAN_2D.MultiDimFit.mH125.7.root // Processing tdrstyle.cc... // Processing contours2D.cxx... // Attaching file higgsCombinehww_RVRF_SCAN_2D.MultiDimFit.mH125.7.root as _file0... contour2D("RF",60,-1,3, "RV",60,-2,5, 1.0,1.0); %ENDCODE%
![]() <--/twistyPlugin--> Higgs Coupling modelsThe previous physics models were considering only some common signal strength modifiers aff Under the hypothesis that the boson observed is a narrow scalar state with the same structure of couplings as the SM Higgs boson, the Higgs production and decay modes decouple and we can parameterize the expected yields in terms of production cross sections and decay widths:where the total width is given by the sum of the partial widths in all SM decay modes, plus possibly some BSM decays The κV, κf benchmark model.In this model, we introduce just two parameters κV, κf that scale the couplings between the Higgs boson and the electroweak gauge bosons and between the Higgs boson and the fermions. The loop-induced effective couplings of the Higgs boson to gluons and to photons are computed in terms of κV, κf by assuming that only the SM particles participate to the loop.
where we have neglected the γγ and Zγ decays that contribute only to the few per-mille level. We can then dividing by the SM width, so that we can express a scale factor for the total width κH in terms of SM BRs Since in the SM for a 125.7 GeV Higgs we have BR(H→VV) = 0.254, to a good approximation Exercises for the reader
<--/twistyPlugin twikiMakeVisibleInline-->
<--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline-->
<--/twistyPlugin--> Implementation in combineThe PhysicsModel that implements the κV and κf model is CvCfHiggs![]() ![]() ![]() ![]() ![]() ![]() Working example: H→WWLet's now investigate the constraints in the (κV,κf) provided by the H→WW analysis separately for 0 jets and 2 jets, and for the combination of the two We can create the three workspaces:text2workspace.py hwwof_0j_shape_8TeV.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww0j.root --for-fits text2workspace.py hwwof_2j_shape_8TeV.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww2j.root --for-fits text2workspace.py comb_hww.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_hww.root --for-fitsand then we run the scans on the Asimov dataset, since they are easier to interpret for X in hww hww0j hww2j; do ( combine -M MultiDimFit cVcF_${X}.root -m 125.7 -n ${X}_CVCF_SCAN_2D_ASIMOV --algo=grid3x3 --points=225 --setPhysicsModelParameterRanges CF=0,2:CV=0,2 -t -1 --expectSignal=1 > log.cvf$X 2>&1 & ); doneThe resulting plot is the following: ![]() <--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_hww.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_hww() { TCanvas *c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_hww.root", "RECREATE"); const int states = 3; const char *state[states] = { "hww", "hww0j", "hww2j" }; for (int i = 0; i < states; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D_ASIMOV.MultiDimFit.mH125.7.root", state[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",state[i])); gFile->Close(); } fOut->cd(); TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS");
TList c68_comb = (TList) fOut->Get("kvkf_hww_c68"); styleMultiGraph(c68_comb, /*color=*/1, /*width=*/3, /*style=*/1); fillMultiGraph(c68_comb, /*color=*/kGray, /*style=*/1001);
TList c68_0j = (TList) fOut->Get("kvkf_hww0j_c68"); styleMultiGraph(c68_0j, /*color=*/4, /*width=*/3, /*style=*/1); fillMultiGraph(c68_0j, /*color=*/4, /*style=*/3005);
TList c68_2j = (TList) fOut->Get("kvkf_hww2j_c68"); fillMultiGraph(c68_2j, /*color=*/2, /*style=*/3004); styleMultiGraph(c68_2j, /*color=*/2, /*width=*/3, /*style=*/1);
c68_2j->Draw("F SAME"); c68_0j->Draw("F SAME"); c68_comb->Draw("F SAME");
c68_2j->Draw("L SAME"); c68_0j->Draw("L SAME"); c68_comb->Draw("L SAME");
TGraph *best = fOut->Get("kvkf_hww_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME");
TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); TLegend *leg = new TLegend(0.20,0.80,0.45,0.95); leg->SetTextFont(42); leg->SetShadowColor(0); leg->SetFillColor(0); leg->AddEntry(best, "Comb.", "P"); leg->AddEntry(c68_comb->At(0), "(68% CL)", "LF"); leg->AddEntry(c68_0j->At(0), "0 jets", "LF"); leg->AddEntry(c68_2j->At(0), "2 jets", "LF"); leg->Draw(); c1->Print("output_kvkf_hww.png"); fOut->Close(); }
%ENDCODE%
<--/twistyPlugin-->The shape of the H→WW + 2 jets can be explained considering that the expected signal is mostly from VBF. As computed in the previous questions from the reader, the expected scale factor for VBF H→WW is as long as κV ≤ 2κf, the κV at the denominator can be neglected, and so the measurement of μ is a measurement of C6: a general six-parameter model for the couplings of a SM-like Higgs bosonThe κV, κf model is quite restrictive, assuming e.g. the same coupling strength modifier for both quarks and leptons, and for both up-type and down-type fermions, assumptions which are not true in many BSM scenarios (e.g. in supersymmetry the couplings to the top and bottom quarks have different modifiers). We will consider therefore a more generic model which relies on a much smaller set of assumptions, namely:
Exercise for the reader
<--/twistyPlugin twikiMakeVisibleInline--> Γtot/ΓtotSM = 0.566¡¤κb2 + 0.2541¡¤κV2 + 0.0851¡¤κg2 + 0.0623¡¤κτ2 + 0.0285¡¤κt2 + 0.00388¡¤κγ2
<--/twistyPlugin--> <--
Implementation in combineThe PhysicsModel that implements the c6 model is C6![]() = 'doHZg': self.doHZg = True def doParametersOfInterest(self): """Create POI out of signal strength and MH""" self.modelBuilder.doVar("kV[1,0.0,2.0]") self.modelBuilder.doVar("ktau[1,0.0,2.0]") self.modelBuilder.doVar("ktop[1,0.0,2.0]") self.modelBuilder.doVar("kbottom[1,0.0,2.0]") self.modelBuilder.doVar("kgluon[1,0.0,2.0]") self.modelBuilder.doVar("kgamma[1,0.0,2.0]") pois = 'kV,ktau,ktop,kbottom,kgluon,kgamma' if self.doHZg: self.modelBuilder.doVar("kZgamma[1,0.0,30.0]") pois + ",kZgamma" if self.floatMass: # [...] else: # [...] self.modelBuilder.doSet("POI",pois)
## Create the SMHiggsBuilder, to get the SM BRs self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup()
def setup(self):
# SM BRs, needed to compute the modified total width as funciton of the partial widths for d in [ "htt", "hbb", "hcc", "hww", "hzz", "hgluglu", "htoptop", "hgg", "hzg", "hmm", "hss" ]: self.SMH.makeBR(d)
## Partial decay witdhs, normalized to the SM one self.modelBuilder.factory_('expr::c6_Gscal_Vectors("@0*@0 * (@1![]() ![]() ![]() ![]() ![]() ![]() = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production = "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay == "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() Working example 1: κγWe will take the full combination but excluding the H→inv search, and extract the κγ modifier, profiling the others. However, to profile them correctly we need to select ranges that are large enough so that the parameter doesn't end on the boundary. We can look at the past CMS Higgs combination, HIG-13-005, to see what was the sensitivity to the other parameters, e.g from the summary plot![]() In addition to the observed outcome, we will compute also the Asimov one text2workspace.py comb.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:c6 -o c6_comb.root --for-fit RANGES="--setPhysicsModelParameterRanges kbottom=0,5:ktop=0,8:ktau=0,4:kgluon=0,4:kgamma=0,3" combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma /// note: 40 points take about 4 minutes (on lxplus) combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV -t -1 --expectSignal=1 /// note: 40 points take about 20 minutes (on lxplus) ![]() Speeding up things: running scans in parallelWhen the model becomes complex, evaluating all the points of the scan in a single job can be inefficient. For this purpose, MultiDimFit allows one to fraction the task by specifying a range of points to run: different ranges can be run in parallel, and at the end the root files can be merged withhadd to get the final one. For example, for the Asimov scan of κγ we could do the 40 points in 5 jobs, doing
combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.0 -t -1 --expectSignal=1 --firstPoint 0 --lastPoint 7 combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.1 -t -1 --expectSignal=1 --firstPoint 8 --lastPoint 15 combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.2 -t -1 --expectSignal=1 --firstPoint 16 --lastPoint 23 combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.3 -t -1 --expectSignal=1 --firstPoint 24 --lastPoint 31 combine c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV.4 -t -1 --expectSignal=1 --firstPoint 32 --lastPoint 39 hadd higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.MultiDimFit.mH125.7.root higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.*.MultiDimFit.mH125.7.rootwhere the first 5 commands can be executed in parallel. For convenience, there is even already a script that takes care of running all the jobs in parallel splitting the points, suitable for running on a multi-core machine: test/parallelScan.py ![]() cp -s ${CMSSW_BASE}/src/HiggsAnalysis/CombinedLimit/test/parallelScan.py . parallelScan.py -j 5 c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgamma --floatOtherPOI=1 $RANGES --points=40 -n C6_comb_SCAN_1D_kgamma_ASIMOV -t -1 --expectSignal=1 hadd higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.MultiDimFit.mH125.7.root higgsCombineC6_comb_SCAN_1D_kgamma_ASIMOV.*.MultiDimFit.mH125.7.rootThe option -j N is used to specify the number of jobs; if not present, it will use one job per CPU of the machine, or half that number if it detects that is being run on a lxplus or cmslpc node (one should be respectful of other users logged in the same machine)
Working example 2: κg vs κb degeneracyThe channels with the best sensitivity to the signal strength are mostly driven by the gluon fusion production cross section (e.g. H→WW and H→ZZ, H→γγboosted H→ττ). Conversely, the H→bb is the dominant decay mode, so that the total width is mostly sensitive to κb. The combination of these two features results in a strong correlation between the uncertainties of the two measurements of κg and κb: a simultaneous increase of both couplings would produce an increase in the gluon fusion production cross section but also a decrease of the branching ratio in the other modes (WW, ZZ, γγ and ττ), leaving the expected σ¡¤BR unchanged, and so provides still a good fit to the data. We can see this feature with a 2D likelihood scan:./parallelScan.py c6_comb.root -m 125.7 -M MultiDimFit --algo=grid -P kgluon -P kbottom --floatOtherPOI=1 $RANGES --points=400 -n C6_comb_SCAN_2D_kgluon_kbottom ![]() The contour in these two parameters is very elongated, exemplifying this correlation. Coupling fits from the CMS combination
<--/twistyPlugin twikiMakeVisibleInline-->
The commands to run to create the workspaces and do the scans are below. Running a 2D scan takes about 2 minutes for a single channel, and about 20 minutes for the full combination
text2workspace.py comb.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_comb.root --for-fit ; for X in h{bb,tt,gg,ww,zz}; do text2workspace.py comb_${X}.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.HiggsCouplings:cVcF -o cVcF_${X}.root --for-fit ; done for X in h{bb,tt,gg,ww,zz,comb}; do combine -M MultiDimFit cVcF_${X}.root -m 125.7 -n ${X}_CVCF_SCAN_2D --algo=grid3x3 --points=225 --setPhysicsModelParameterRanges CV=0,2:CF=0,2; doneA macro to make the summary plot is below: %CODE{"cpp"}% // run with root.exe -b -l tdrstyle.cc contours2D.cxx plotcvcf_all.cxx -q void fillMultiGraph(TList tmg, int fillColor, int fillStyle) { for (int i = 0; i < tmg->GetSize(); ++i) { TGraph *g = (TGraph) tmg->At(i); g->SetFillColor(fillColor); g->SetFillStyle(fillStyle); } } void plotcvcf_all() { TCanvas c1 = new TCanvas("c1","c1"); TFile *fOut = TFile::Open("kvkf_plots_all.root", "RECREATE"); const int decays = 6; const char *decay[decays] = { "hbb", "htt", "hzz", "hww", "hgg", "comb" }; const int color[decays] = { 67 , 221 , 2 , 4 , 209 , 1 }; const int fill[decays] = { 228 , 223 , 208 , 216 , 211 , 203 }; TList *c68[decays] = { 0,0,0,0,0, 0 }; for (int i = 0; i < decays; ++i) { gFile = TFile::Open(Form("higgsCombine%s_CVCF_SCAN_2D.MultiDimFit.mH125.7.root", decay[i])); contour2D("CV",45,0.,2., "CF",45,0.,2., 1.0,1.0, fOut,Form("kvkf_%s",decay[i])); gFile->Close(); } fOut->cd(); fOut->ls(); for (int i = 0; i < decays; ++i) { c68[i] = (TList) fOut->Get(Form("kvkf_%s_c68",decay[i])); styleMultiGraph(c68[i], color[i], 3, 1); fillMultiGraph(c68[i], fill[i], 1001); } TH2F *frame = new TH2F("frame",";#kappa_{V};#kappa_{f}", 45,0.,2., 45,0.,2.); frame->GetXaxis()->SetDecimals(1); frame->GetYaxis()->SetDecimals(1); frame->Draw("AXIS"); for (int i = 0; i < decays; ++i) { c68[i]->Draw("F SAME"); } for (int i = 0; i < decays; ++i) { c68[i]->Draw("L SAME"); } TGraph *best = fOut->Get("kvkf_comb_best"); best->SetMarkerColor(1); best->SetMarkerStyle(34); // apparently ROOT forgets this best->SetMarkerSize(2.0); // so we do this again best->Draw("P SAME"); TMarker m; m.SetMarkerSize(3.0); m.SetMarkerColor(97); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); m.SetMarkerSize(1.8); m.SetMarkerColor(89); m.SetMarkerStyle(33); m.DrawMarker(1.0,1.0); c1->Print("output_kvkf_all.png"); fOut->Close(); } %ENDCODE% ![]() <--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline--> combine c6_comb.root -M MultiDimFit -m 125.7 $RANGES >>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods. Set Range of Parameter kbottom To : (0,5) Set Range of Parameter ktop To : (0,8) Set Range of Parameter ktau To : (0,4) Set Range of Parameter kgluon To : (0,4) Set Range of Parameter kgamma To : (0,3) Computing limit starting from observation --- MultiDimFit --- best fit parameter values: kV : +1.044 ktau : +1.274 ktop : +2.053 kbottom : +1.401 kgluon : +1.039 kgamma : +1.143 Done in 0.11 min (cpu), 0.11 min (real) combine c6_comb.root -M MultiDimFit --algo=singles -m 125.7 $RANGES >>> including systematics >>> method used to compute upper limit is MultiDimFit >>> random number generator seed is 123456 ModelConfig 'ModelConfig' defines more than one parameter of interest. This is not supported in some statistical methods. Set Range of Parameter kbottom To : (0,5) Set Range of Parameter ktop To : (0,8) Set Range of Parameter ktau To : (0,4) Set Range of Parameter kgluon To : (0,4) Set Range of Parameter kgamma To : (0,3) Computing limit starting from observation Running Minos for POI Real time 0:00:11, CP time 11.700 Running Minos for POI Real time 0:00:11, CP time 11.700 Running Minos for POI Real time 0:00:11, CP time 11.700 Running Minos for POI Real time 0:00:11, CP time 11.700 Running Minos for POI Real time 0:00:11, CP time 11.700 Running Minos for POI Real time 0:00:11, CP time 11.700 --- MultiDimFit --- best fit parameter values and profile-likelihood uncertainties: kV : +1.044 -0.234/+0.204 (68%) ktau : +1.274 -0.497/+0.504 (68%) ktop : +2.053 -0.897/+1.702 (68%) kbottom : +1.401 -0.674/+1.227 (68%) kgluon : +1.039 -0.324/+0.811 (68%) kgamma : +1.143 -0.339/+0.388 (68%) Done in 1.55 min (cpu), 1.64 min (real) <--/twistyPlugin--> Going beyond: introducing H→invis decays into C6The C6 model assumes that there are no BSM decays of the Higgs. We will now take one step forward, and allow the possibility of Higgs decays into invisible particles, relying on the searches in the VBF mode to constrain this. We will introduce a new parameter BRInv, representing the BR of the Higgs boson decays to invisible particles.
<--/twistyPlugin twikiMakeVisibleInline--> The total width is given by the sum of the partial widths of the SM-like decays (dependent on κ) plus the invisible one, and by definition and consequently dividing by the SM total width we can express this in terms of SM BRs <--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"python"}% class C6I(SMLikeHiggsModel): "as C6 but including direct searches for H->invis" def __init__(self): SMLikeHiggsModel.__init__(self) # not using 'super(x,self).__init__' since I don't understand it self.floatMass = False self.doHZg = False def setPhysicsOptions(self,physOptions): for po in physOptions: if po.startswith("higgsMassRange="): self.floatMass = True self.mHRange = po.replace("higgsMassRange=","").split(",") print 'The Higgs mass range:', self.mHRange if len(self.mHRange) = 2: raise RuntimeError, "Higgs mass range definition requires two extrema" elif float(self.mHRange[0]) >= float(self.mHRange[1]): raise RuntimeError, "Extrama for Higgs mass range defined with inverterd order. Second must be larger the first" if po = 'doHZg': self.doHZg = True def doParametersOfInterest(self): self.modelBuilder.doVar("kV[1,0.0,2.0]") self.modelBuilder.doVar("ktau[1,0.0,2.0]") self.modelBuilder.doVar("ktop[1,0.0,2.0]") self.modelBuilder.doVar("kbottom[1,0.0,2.0]") self.modelBuilder.doVar("kgluon[1,0.0,2.0]") self.modelBuilder.doVar("kgamma[1,0.0,2.0]") self.modelBuilder.doVar("BRInv[0,0,1]") pois = 'kV,ktau,ktop,kbottom,kgluon,kgamma,BRInv' if self.doHZg: self.modelBuilder.doVar("kZgamma[1,0.0,30.0]") pois + ",kZgamma" if self.floatMass: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setRange(float(self.mHRange[0]),float(self.mHRange[1])) self.modelBuilder.out.var("MH").setConstant(False) else: self.modelBuilder.doVar("MH[%s,%s]" % (self.mHRange[0],self.mHRange[1])) self.modelBuilder.doSet("POI",pois+',MH') else: if self.modelBuilder.out.var("MH"): self.modelBuilder.out.var("MH").setVal(self.options.mass) self.modelBuilder.out.var("MH").setConstant(True) else: self.modelBuilder.doVar("MH[%g]" % self.options.mass) self.modelBuilder.doSet("POI",pois) self.SMH = SMHiggsBuilder(self.modelBuilder) self.setup()
def setup(self):
# SM BR for d in [ "htt", "hbb", "hcc", "hww", "hzz", "hgluglu", "htoptop", "hgg", "hzg", "hmm", "hss" ]: self.SMH.makeBR(d)
## total witdhs, normalized to the SM one self.modelBuilder.factory_('expr::c6i_Gscal_Vectors("@0*@0 * (@1![]() ![]() ![]() ![]() ![]() ![]() = None: XSscal = "kgluon" if production in ["WH","ZH","VH","qqH"]: XSscal = "kV" if production = "ttH": XSscal = "ktop" BRscal = "hgg" if decay in ["hbb", "htt", "hinv"]: BRscal = decay if decay in ["hww", "hzz"]: BRscal = "hvv" if self.doHZg and decay == "hzg": BRscal = "hzg" self.modelBuilder.factory_('expr::%s("@0*@0 * @1![]() <--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline--> text2workspace.py hinv_vbf_8TeV.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_vbfinv.root --for-fit combine -m 125.7 -M MultiDimFit c6i_vbfinv.root -n C6I_vbfinv_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=400 --algo=grid3x3 $RANGES text2workspace.py comb_hww.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_hww.root --for-fit ./parallelScan.py -M MultiDimFit c6i_hww.root -n C6I_hww_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=100 --algo=grid3x3 $RANGES hadd -f higgsCombineC6I_hww_SCAN_2D_kV_BRInv.MultiDimFit.mH125.7.root higgsCombineC6I_hww_SCAN_2D_kV_BRInv.*.MultiDimFit.mH125.7.root HWW="hww0j_8TeV=hwwof_0j_shape_8TeV.txt hww2j_8TeV=hwwof_2j_shape_8TeV.txt" HIN="hinv_8TeV=hinv_vbf_8TeV.txt" combineCards.py -S $HWW $HIN > comb_test2.txt text2workspace.py comb_test2.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_test2.root --for-fit ./parallelScan.py -m 125.7 -M MultiDimFit c6i_test2.root -n C6I_test2_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=0 --points=100 --algo=grid3x3 $RANGES hadd -f higgsCombineC6I_test2_SCAN_2D_kV_BRInv.MultiDimFit.mH125.7.root higgsCombineC6I_test2_SCAN_2D_kV_BRInv.*.MultiDimFit.mH125.7.root
<--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline--> text2workspace.py comb.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_comb.root --for-fit; text2workspace.py combx.txt -m 125.7 -P HiggsAnalysis.CombinedLimit.LOFullParametrization:c6i -o c6i_combx.root --for-fit; ./parallelScan.py -m 125.7 -M MultiDimFit c6i_comb.root -n C6I_comb_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=1 --points=400 --algo=grid $RANGES --hadd ./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_2D_kV_BRInv -P kV -P BRInv --floatOtherPOI=1 --points=400 --algo=grid $RANGES ---hadd ![]() The solid-filled areas are the 68% and 95% CL regions for the combination not including the direct search for H→invis, the hatched areas are the 68% and 95% CL regions including it. You can appreciate how without the constraint there is a full degeneracy between the measured coupling and the invisible BR, which is then broken by the direct search. <--/twistyPlugin-->
<--/twistyPlugin twikiMakeVisibleInline-->
Best fit and one-dimensional uncertainties (observed): takes O(1h)
combine -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_comb_SINGLES --algo=singles $RANGES --robustFit=1 --- MultiDimFit --- best fit parameter values and profile-likelihood uncertainties: kV : +1.146 -0.207/+0.199 (68%) ktau : +1.399 -0.466/+0.428 (68%) ktop : +2.000 -0.683/+0.000 (68%) kbottom : +1.428 -0.537/+0.572 (68%) kgluon : +1.075 -0.265/+0.407 (68%) kgamma : +1.231 -0.306/+0.345 (68%) BRInv : +0.177 -0.177/+0.157 (68%) Done in 54.14 min (cpu), 55.02 min (real)One-dimensional expected uncertainties: one parameter at a time in parallel, since it takes O(1h) per parameter. for D in k{V,gamma,gluon,bottom,top,tau} BRInv; do ( combine -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_comb_SINGLES_${D}_Asimov --algo=singles $RANGES --robustFit=1 -t -1 -P $D --floatOtherPOI=1 --expectSignal=1 2>&1 | tee log.C6I_Asimov_$D > /dev/null & ); done tail -n 2 log.C6I_Asimov_* ==> log.C6I_Asimov_BRInv <== BRInv : +0.039 -0.039/+0.190 (68%) Done in 58.08 min (cpu), 65.68 min (real) ==> log.C6I_Asimov_kV <== kV : +1.020 -0.201/+0.196 (68%) Done in 62.85 min (cpu), 70.62 min (real) ==> log.C6I_Asimov_kbottom <== kbottom : +1.020 -0.433/+0.585 (68%) Done in 66.84 min (cpu), 74.38 min (real) ==> log.C6I_Asimov_kgamma <== kgamma : +1.021 -0.225/+0.255 (68%) Done in 61.08 min (cpu), 68.59 min (real) ==> log.C6I_Asimov_kgluon <== kgluon : +1.017 -0.241/+0.457 (68%) Done in 74.38 min (cpu), 81.99 min (real) ==> log.C6I_Asimov_ktau <== ktau : +1.022 -0.366/+0.323 (68%) Done in 74.65 min (cpu), 82.36 min (real) ==> log.C6I_Asimov_ktop <== ktop : +1.019 -1.019/+0.750 (68%) Done in 57.58 min (cpu), 65.44 min (real)One-dimensional likelihood scans (observed and expected) for D in k{V,gamma,gluon,bottom,top,tau} BRInv; do ./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_1D_$D -P ${D} --floatOtherPOI=1 $RANGES --algo=grid --points 50 --hadd; ./parallelScan.py -m 125.7 -M MultiDimFit c6i_combx.root -n C6I_combx_SCAN_1D_${D}_ASIMOV -P ${D} --floatOtherPOI=1 $RANGES --algo=grid --points 30 -t -1 --expectSignal=1 --hadd; doneResults:
<--/twistyPlugin twikiMakeVisibleInline-->
%CODE{"cpp"}% void plot1ds_c6i() { TCanvas c1 = new TCanvas("c1","c1"); gStyle->SetOptStat(0); TFile *fOut = TFile::Open("plot1d_c6i.root", "RECREATE"); const int params = 7; const char *param[params] = { "kV", "kbottom", "ktop", "ktau", "kgluon", "kgamma", "BRInv" }; const float min[params] = { 0, 0, 0, 0, 0, 0, 0, }; const float max[params] = { 2, 5, 8, 4, 4, 3, 1, }; const char *label[params] = { "#kappa_{V}", "#kappa_{b}", "#kappa_{t}", "#kappa_{#tau}", "#kappa_{g}", "#kappa_{#gamma}","BR_{Inv}" }; for (int i = 0; i < params; ++i) { cout << param[i] << endl; gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s.MultiDimFit.mH125.7.root", param[i])); TTree *limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph gO = gROOT->FindObject("Graph")->Clone(); gO->Sort(); gFile->Close(); gFile = TFile::Open(Form("higgsCombineC6I_combx_SCAN_1D_%s_ASIMOV.MultiDimFit.mH125.7.root", param[i])); limit = (TTree) gFile->Get("limit"); limit->Draw(Form("2*deltaNLL:%s",param[i]),""); TGraph *gA = gROOT->FindObject("Graph")->Clone(); gA->Sort(); gFile->Close();
gO->SetName(Form("c6i_1d_%s_obs", param[i])); gA->SetName(Form("c6i_1d_%s_asimov",param[i])); fOut->WriteTObject(gO); fOut->WriteTObject(gA);
TH1F *frame = new TH1F("frame",Form(";%s;- 2 #Delta ln L",label[i]), 10, min[i], max[i]); frame->Draw(); frame->GetYaxis()->SetRangeUser(0,7); gO->SetLineWidth(3); gO->SetLineColor(1); gA->SetLineWidth(2); gA->SetLineColor(1); gA->SetLineStyle(2); gO->Draw("C SAME"); gA->Draw("C SAME");
c1->Print(Form("output_C6I_1D_%s.png", param[i])); delete frame; } fOut->Close(); }
%ENDCODE%
<--/twistyPlugin--> <--/twistyPlugin--> Higgs JCPIntroductionThe narrow resonance with an invariant mass of approximately 125 GeV that has been discovered by ATLAS and CMS, are mainly via diboson decay modes γγ, WW and ZZ. The determination of the quantum numbers of this Higgs-like particle, such as its spin and parity, is currently one of the most important tasks in high energy physics and thus is important for future studies.KinematicsThe scattering amplitudes of Higgs or any exotic boson decays to the WW final state have been performed in several literatures. It was demonstrated that the WW chaneel has a good sensitivity to distinguish between the standard model Higgs boson hypothesis and a spin-2 resonance, which couples to the dibosons through minimal couplings, referred to as 2+. Some sensitivity in this final state can also distinguish between the min0+ and the pseudoscalar 0− boson hypotheses. The typical experimental observables are m(ll), mT(ll,ETmiss), and Δφ(ll). The tdefinition of transverse mass mT(ll,ETmiss) used in this analysis ismT(ll,ETmiss) := sqrt( 2 * pT(ll) * ETmiss * (1 - cos Δ#phi;(ll, ETmiss) )
Exercise for the readerTwo datasets of simulated X→WW→2l2ν events have been produced, with X being a SM Higgs and spin 2 graviton, and a very basic event selection has already been applied to retain only the events with at least two high-pt leptons, and no hadronic jets. What you can do is to run on those simulated events and make plots of the kinematic variables like the ones above (or even others, e.g. ΔR(ll)). You can do this exercise either starting from simple ROOT trees, or from PAT files (for the latter, use CMSSW_5_3_14), using any of the analysis tools you should have learned from the pre-exercises (root, fwlite, python, ...).
cmsStageIn /store/cmst3/group/das2014/HiggsProperties/FILE.root . (or, temporarily, copy them from /afs/cern.ch/work/g/gpetrucc/GENS/CMSSW_5_3_14/src/ ). The files have been produced from these two datasets:
/afs/cern.ch/user/g/gpetrucc/public/CMSDAS-2014-CERN-toy_HWW.py (note: this is just a toy example, the real H→WW selection has many more selection requirements, and use more objects and tools... see HIG-13-023 for more information)
<--/twistyPlugin twikiMakeVisibleInline-->
Here are many possible solutions, depending on whether you run on the Trees or on the PAT, and what software you want to use. Here we will give four.
Using Trees, in python %CODE{"python"}% from math import * import sys, ROOT from DataFormats.FWLite import Handle, Events def fillPlots(tree,histos): l1 = ROOT.TLorentzVector() l2 = ROOT.TLorentzVector() met = ROOT.TLorentzVector() print "Making plots for one sample: ", for i in xrange(tree.GetEntries()): # read tree.GetEntry(i) l1.SetPtEtaPhiM(tree.lep1_pt, tree.lep1_eta, tree.lep1_phi, tree.lep1_mass) l2.SetPtEtaPhiM(tree.lep2_pt, tree.lep2_eta, tree.lep2_phi, tree.lep2_mass) met.SetPtEtaPhiM(tree.met_pt, 0.0, tree.met_phi, 0.0) # -- m(ll) p4ll = l1 + l2 histos['mll'].Fill(p4ll.M()) # -- delta phi(ll) dphi = abs(l1.Phi() - l2.Phi()) if dphi > pi: dphi = 2*pi-dphi histos['dphill'].Fill(dphi) # -- mT mT = sqrt(2*p4ll.Pt()*met.Pt()*(1-cos(p4ll.Phi()-met.Phi()))) histos['mT'].Fill(mT) # -- 2d histos['m2d'].Fill(mT,p4ll.M()) if (i+1) % 1000 == 0: sys.stdout.write("."); sys.stdout.flush() print " done." def stackPlots(plots): c1 = ROOT.TCanvas("c1") for var,plotsm in plots["sm"].iteritems(): plot2p = plots["2p"][var] if "TH1F" in plotsm.ClassName(): # common things for p in plotsm, plot2p: p.SetLineWidth(2) p.Scale(1.0/p.Integral()) # different things plotsm.SetLineColor(1) plot2p.SetLineColor(2) # plot plotsm.Draw("HIST ") plot2p.Draw("HIST SAME") plotsm.GetYaxis().SetRangeUser(0, 1.4*max(plotsm.GetMaximum(),plot2p.GetMaximum())) # legend leg = ROOT.TLegend(0.2, 0.77, 0.55, 0.92) leg.SetFillColor(0) leg.SetTextFont(42) leg.SetTextSize(0.05) leg.AddEntry(plotsm, "SM Higgs", "L") leg.AddEntry(plot2p, "2^{+}_{m} Graviton", "L") leg.Draw(); # print c1.Print("plot_jcp_"+var+".png") elif "TH2F" in plotsm.ClassName(): plotsm.Draw("COLZ") c1.Print("plot_jcp_"+var+"_sm.png") plot2p.Draw("COLZ") c1.Print("plot_jcp_"+var+"_2p.png") ROOT.gROOT.ProcessLine(".x tdrstyle.cc") ROOT.gROOT.SetBatch(1) ROOT.gStyle.SetOptStat(0) output = ROOT.TFile("plots.root","RECREATE") plots = {} for spin in "2p", "sm": plots[spin] = { "mll" :ROOT.TH1F("mll_" +spin, ";m(ll) (GeV)",40,0,120), "mT" :ROOT.TH1F("mT_" +spin, ";m_{T}(ll,E_{T}^{miss}) (GeV)",40,0,200), "dphill":ROOT.TH1F("dphill_"+spin, ";#Delta#phi(ll) (rad)",31,0,pi), "m2d" :ROOT.TH2F("m2d_"+spin, ";m_{T}(ll,E_{T}^{miss}) (GeV);m(ll) (GeV)",6,60,140,5,12,100), } smhiggs_file = ROOT.TFile.Open("HWW_Tree_SMHiggs.root") graviton_file = ROOT.TFile.Open("HWW_Tree_Graviton.root") smhiggs = smhiggs_file.Get("tree/probe_tree") graviton = graviton_file.Get("tree/probe_tree") fillPlots(graviton, plots["2p"]) fillPlots(smhiggs, plots["sm"]) stackPlots(plots) %ENDCODE% Using PAT Objects, in python (with FWLite) %CODE{"python"}% from math import * import sys, ROOT from DataFormats.FWLite import Handle, Events def fillPlots(events,histos): muons = Handle("std::vector") electrons = Handle("std::vector") mets = Handle("std::vector") print "Making plots for one sample: ", for i,event in enumerate(events): event.getByLabel("cleanPatMuons", muons) event.getByLabel("cleanPatElectrons", electrons) event.getByLabel("patMETs", mets) # list of leptons sorted by pT leptons = [ mu.p4() for mu in muons.product() ] + [ el.p4(el.candidateP4Kind()) for el in electrons.product() ] leptons.sort(key = lambda lepton : lepton.Pt(), reverse=True) # missing energy met = mets.product()[0].p4() # reconstruct the various kinematic variables # -- m(ll) p4ll = leptons[0] + leptons[1] histos['mll'].Fill(p4ll.mass()) # -- delta phi(ll) dphi = abs(leptons[0].Phi() - leptons[1].Phi()) if dphi > pi: dphi = 2*pi-dphi histos['dphill'].Fill(dphi) # -- mT mT = sqrt(2*p4ll.Pt()*met.Pt()*(1-cos(p4ll.Phi()-met.Phi()))) histos['mT'].Fill(mT) # -- 2d histos['m2d'].Fill(mT,p4ll.M()) if (i+1) % 1000 == 0: sys.stdout.write("."); sys.stdout.flush() print " done." def stackPlots(plots): c1 = ROOT.TCanvas("c1") for var,plotsm in plots["sm"].iteritems(): plot2p = plots["2p"][var] if "TH1F" in plotsm.ClassName(): # common things for p in plotsm, plot2p: p.SetLineWidth(2) p.Scale(1.0/p.Integral()) # different things plotsm.SetLineColor(1) plot2p.SetLineColor(2) # plot plotsm.Draw("HIST ") plot2p.Draw("HIST SAME") plotsm.GetYaxis().SetRangeUser(0, 1.4*max(plotsm.GetMaximum(),plot2p.GetMaximum())) # legend leg = ROOT.TLegend(0.2, 0.77, 0.55, 0.92) leg.SetFillColor(0) leg.SetTextFont(42) leg.SetTextSize(0.05) leg.AddEntry(plotsm, "SM Higgs", "L") leg.AddEntry(plot2p, "2^{+}_{m} Graviton", "L") leg.Draw(); # print c1.Print("plot_jcp_"+var+".png") elif "TH2F" in plotsm.ClassName(): plotsm.Draw("COLZ") c1.Print("plot_jcp_"+var+"_sm.png") plot2p.Draw("COLZ") c1.Print("plot_jcp_"+var+"_2p.png") ROOT.gROOT.ProcessLine(".x tdrstyle.cc") ROOT.gROOT.SetBatch(1) ROOT.gStyle.SetOptStat(0) output = ROOT.TFile("plots.root","RECREATE") plots = {} for spin in "2p", "sm": plots[spin] = { "mll" :ROOT.TH1F("mll_" +spin, ";m(ll) (GeV)",40,0,120), "mT" :ROOT.TH1F("mT_" +spin, ";m_{T}(ll,E_{T}^{miss}) (GeV)",40,0,200), "dphill":ROOT.TH1F("dphill_"+spin, ";#Delta#phi(ll)",31,0,pi), "m2d" :ROOT.TH2F("m2d_"+spin, ";m_{T}(ll,E_{T}^{miss}) (GeV);m(ll) (GeV)",6,60,140,5,12,100), } smhiggs = Events("HWW_SMHiggs.root") graviton = Events("HWW_Graviton.root") fillPlots(graviton, plots["2p"]) fillPlots(smhiggs, plots["sm"]) stackPlots(plots) %ENDCODE% Using Trees, in C++ (with ROOT) %CODE{"cpp"}% #include #include #include #include #include #include #include #include #include #include For this to work you must always compile the macro (i.e. run root.exe -b -l -q makePlots_PAT.cxx++ ), and you need to have these lines, or something equivalent, in your rootlogon
%CODE{"cpp"}% //// == rootlogon === { TString makeshared(gSystem->GetMakeSharedLib()); TString dummy = makeshared.ReplaceAll("-W ", ""); dummy.ReplaceAll("g++ ","g++ -std=c++0x "); dummy.ReplaceAll("-Wshadow "," "); dummy.ReplaceAll("-Wall ","-Wall -Wno-deprecated-declarations "); gSystem->SetMakeSharedLib(dummy); gSystem->Load("libGenVector.so"); // strangely this is needed if (getenv("CMSSW_BASE")) { gSystem->Load("libFWCoreFWLite.so"); gSystem->Load("libDataFormatsFWLite.so"); gSystem->SetIncludePath("-I$ROOFITSYS/include"); AutoLibraryLoader::enable(); } }
%ENDCODE%
%CODE{"cpp"}% //// = makePlots_PAT.cxx++ = #include #include #include #include #include #include #include #include #include #include #include <--/twistyPlugin--> StatisticsThe statistics of hypotheses testing is based on the likelihood ratio between the hypotheses being tested.Likelihood ratioConsider the null (N) and alternate (A) hypotheses which correspond in this case to the SMH-like 0++background (N) and the alternative JCP+background (A). The likelihood ratio is defined as q = -2 ln( L(x|A) / L(x|N) ), where x represents a set of events. This set of events, x, can be:
![]() ![]() Measures of separation and the CLs criteriumGiven the distributions above, one can then characterise the expectations and observation in a number of ways:
Pre-fit and post-fit expectationsA final complication relates to the normalisation of the processes entering the N and A hypotheses, especially the signal normalisation. We consider the "pre-fit" and "post-fit" cases:
![]() ![]() ![]() ![]() ![]() Exercise for the reader
What does it mean to show -0.9σ in the Table above?Hint: how do you describe the case in which P( q < qobs | N ) < 0.5? <--/twistyPlugin twikiMakeVisibleInline--> The answer is deceivingly simple: since P( q < qobs | N ) < 0.5, one needs to have signed deviations to note that qobs is to the left of the median of P( q | N ). In this case, since P( q < qobs | N ) < 0.5 and P( q > qobs | A ) > 0.5, it means that the observation is somewhere in between the qexp(A) and qexp(N) as can be seen in the following figure:
![]() ![]() <--/twistyPlugin--> PracticeFor the spin-2 (2+m in this case) vs spin-0 (0+) hypotheses test, the full H→WW→2a„¡°2ν analysis includes 4 channels corresponding to {0-jet,1-jet}⊗{7 TeV, 8 TeV}. For simplicity we will perform the exercise on the 0-jet 8 TeV data card.head jcp_hwwof_0j_8TeV.txt [...] Observation 5747 shapes * * jcp_hwwof_0j.input_8TeV.root histo_$PROCESS histo_$PROCESS_$SYSTEMATIC shapes data_obs * jcp_hwwof_0j.input_8TeV.root histo_Data bin j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev j0of8tev process qqbarH_ALT ggH_ALT ZH WH qqH ggH qqWW ggWW VV Top Zjets WjetsE Wgamma Wg3l Ztt WjetsM qqWW2j process -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 rate 129.869 146.427 1.728 5.833 3.033 234.884 3973.542 210.835 131.912 499.856 0.000 279.465 115.142 166.923 46.430 327.212 0.288You will note that besides the usual SM signal contributions (ggH, qqH, ZH, and WH) there are two new processes (qqbarH_ALT and ggH_ALT) corresponding to the production of 2+m via quark-antiquark annihilation and gluon-fusion. Let's look at the shape of the templates: root -l jcp_hwwof_0j.input_8TeV.root%CODE{"cpp"}% //Compare main signal and background histo_ggH->GetXaxis()->SetRangeUser(-1,1) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_qqWW->SetLineColor(kGray) histo_qqWW->DrawNormalized("Hsame") %ENDCODE% ![]() ![]() ![]() ![]() root -l jcp_hwwof_0j.input_8TeV.root%CODE{"cpp"}% //Define a region where it is easy to see the differences between the alternate signals histo_ggH->;GetXaxis()->SetRangeUser(-1,0) histo_ggH->SetLineColor(kBlack) histo_ggH->DrawNormalized("H") histo_ggH_ALT->SetLineColor(kBlue) histo_ggH_ALT->DrawNormalized("Hsame") histo_qqbarH_ALT->SetLineColor(kRed) histo_qqbarH_ALT->DrawNormalized("Hsame") %ENDCODE% ![]() ![]() ![]() Building a nested modelThe hypothesis test is performed by constructing a nested model which has two relevant parameters:
![]() text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose -o jcp_hww.root MH (not there before) will be assumed to be 125.7 Process ggH_ALT will get norm x Process ZH will get norm not_x Process qqH will get norm not_x Process qqbarH_ALT will get norm x Process WH will get norm not_x Process ggH will get norm not_x Process ggH_ALT will get norm x Process ZH will get norm not_x Process qqH will get norm not_x Process qqbarH_ALT will get norm x Process WH will get norm not_x Process ggH will get norm not_xSee how the ggH_ALT and qqbarH_ALT are both being scaled with x? That is not what we want, since we would like each of the processes to be differentiated via fqq. To do this properly we need to tell the model to also include fqq: text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose --PO fqqIncluded -o jcp_hww.root Will consider fqq = fraction of qqH in Alt signal (signal strength will be left floating) MH (not there before) will be assumed to be 125.7 Process ggH_ALT will get scaled by r_times_x_times_not_fqq Process ZH will get scaled by r_times_not_x Process qqH will get scaled by r_times_not_x Process qqbarH_ALT will get scaled by r_times_x_times_fqq Process WH will get scaled by r_times_not_x Process ggH will get scaled by r_times_not_x Process ggH_ALT will get scaled by r_times_x_times_not_fqq Process ZH will get scaled by r_times_not_x Process qqH will get scaled by r_times_not_x Process qqbarH_ALT will get scaled by r_times_x_times_fqq Process WH will get scaled by r_times_not_x Process ggH will get scaled by r_times_not_xIn this case you see how ggH_ALT is scaled with r_times_x_times_not_fqq, i.e. r*x*(1-fqq), and qqbarH_ALT is scaled with r_times_x_times_fqq, i.e. r*x*fqq. Running toysLet's dive immediately into throwing toys for the three cases of interest:
GEN_OPTS="-m 125.7 jcp_hww.root --seed 8192" WHAT_OPTS="--testStat=TEV --singlePoint 1 --saveHybridResult" HOW_OPTS="-M HybridNew -T 500 -i 2 --fork 6 --clsAcc 0 --fullBToys" combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --generateExt=1 --generateNuis=0 --expectedFromGrid 0.5 -n jcp_hww_pre-fit_exp &> jcp_hww_pre-fit_exp.log combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist --expectedFromGrid 0.5 -n jcp_hww_post-fit_exp &> jcp_hww_post-fit_exp.log combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist -n jcp_hww_post-fit_obs &> jcp_hww_post-fit_obs.log grep -H "^CLs" jcp_hww_{pre,post}*.logWhile the commands are running (they'll take a few minutes), let's decode the command lines (taken from "combine --help"). There is a first set of common parameters which relate to what we are doing:
jcp_hww_pre-fit_exp.log:CLs = 0.0582329 +/- 0.0108136 jcp_hww_post-fit_exp.log:CLs = 0.206827 +/- 0.0203793 jcp_hww_post-fit_obs.log:CLs = 0.089172 +/- 0.0171273You can see two things:
r at a fixed value of the mixing x corresponding to the two hypotheses
text2workspace.py jcp_hwwof_0j_8TeV.txt -P HiggsAnalysis.CombinedLimit.HiggsJPC:twoHypothesisHiggs -m 125.7 --PO verbose --PO fqqIncluded -o jcp_hww_2poi.root --PO muAsPOI combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P r --floatOtherPOI=0 --setPhysicsModelParameters x=0 -n MU_SCALAR --- MultiDimFit --- best fit parameter values: r : +0.673 combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P r --floatOtherPOI=0 --setPhysicsModelParameters x=1 -n MU_SPIN2 best fit parameter values: r : +0.988This implies that our post-fit expectation for the signal yield is less than the SM Higgs prediction, especially for the spin 0 case, and consequently also our post-fit expectation for the separation between the two hypotheses is worse (with less events, it's harder to separate them). While the approach with toys and CLs is statistically more sound (or at least more in line with the conventions used in particle physics), we can see how well the two models fit the data also by looking at the likelihood as function of x , while profiling r .
combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P x --floatOtherPOI=1 --algo=grid -n JCP_X_SCAN --points=200 combine jcp_hww_2poi.root -M MultiDimFit -m 125.7 -P x --floatOtherPOI=1 --algo=grid -n JCP_X_SCAN_Asimov --points=200 -t -1 --expectSignal=1 --setPhysicsModelParameters x=0,r=1 ![]() From there you can see that the data prefers the scalar hypothesis (x=0), and in particular the events are even more scalar-like than what expected on average for a scalar (the likelihood approaches zero with a non-zero slope). You can also see that for larger x the a-priori expected curve grows faster for increasing x (i.e. it has a smaller that the a-priori sensitivity is better) but that is in part compensated by the fact that the observed starts faster. For a gaussian likelihood, a difference in 2*Δ log L of about 3 would correspond to a p-value of about 10%, which is in the same ballpark as the CLs number; a precise comparison cannot be made, since the two kinds of tests are different and since we're considering a variable at the boundary of the interval in which it is defined (and so it is not really expected to behave as a gaussian) In all the results above, fqq is not floated: it is a constant parameter set at zero. This means that everything we have done up to now tests gg→2+m vs 0+. Exercise for the reader
Redo the same tests for the fqq=1 case?Hint: look at "combine --help" and search for "PhysicsModelParameter". <--/twistyPlugin twikiMakeVisibleInline--> The option allowing to set a model parameter value is --setPhysicsModelParameters.
To re-run CLs, one can do
combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --generateExt=1 --generateNuis=0 --expectedFromGrid 0.5 -n jcp_hww_pre-fit_exp1 --setPhysicsModelParameters fqq=1 &> jcp_hww_pre-fit_exp1.log combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist --expectedFromGrid 0.5 -n jcp_hww_post-fit_exp1 --setPhysicsModelParameters fqq=1 &> jcp_hww_post-fit_exp1.log combine $GEN_OPTS $WHAT_OPTS $HOW_OPTS --frequentist -n jcp_hww_post-fit_obs --setPhysicsModelParameters fqq=1 &> jcp_hww_post-fit_obs1.logThe output is for the two cases is:
![]() <--/twistyPlugin--> PlottingFor simplicity, let's work with the results of the post-fit observed run above which are saved inhiggsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root .
To do this we need to use a macro![]() root -l -q -b higgsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root "${CMSSW_BASE}/src/HiggsAnalysis/CombinedLimit/test/plotting/hypoTestResultTree.cxx(\"jcp_hww_post-fit_obs.qvals.root\",125.7,1,\"x\")" &> jcp_hww_post-fit_obs.qvals.log cat jcp_hww_post-fit_obs.qvals.log Attaching file higgsCombinejcp_hww_post-fit_obs.HybridNew.mH125.7.8192.root as _file0... Processing /Users/adavid/workspace/CMSSW611-master/src/HiggsAnalysis/CombinedLimit/test/plotting/hypoTestResultTree.cxx("jcp_hww_post-fit_obs.qvals.root",125.7,1,"x")... - HypoTestResult_mh125.7_x1_2618602738 Saved test statistics distributions for 996 signal toys and 996 background toys to jcp_hww_post-fit_obs.qvals.root.You can see that the macro found 996 toys in total, which correctly matches the parallel work of 6 processors (166*6=996 while 167*6=1002). Open the output file jcp_hww_post-fit_obs.qvals.root :
root -l jcp_hww_post-fit_obs.qvals.rootand let's explore: %CODE{"cpp"}% q->Print() %ENDCODE% ****************************************************************************** *Tree :q : Test statistics * *Entries : 1993 : Total = 42889 bytes File Size = 8799 * * : : Tree compression factor = 4.95 * ****************************************************************************** *Br 0 :q : q/F * *Entries : 1993 : Total Size= 8495 bytes File Size = 7539 * *Baskets : 1 : Basket Size= 32000 bytes Compression= 1.07 * *............................................................................* *Br 1 :mh : mh/F * *Entries : 1993 : Total Size= 8500 bytes File Size = 143 * *Baskets : 1 : Basket Size= 32000 bytes Compression= 56.21 * *............................................................................* *Br 2 :weight : weight/F * *Entries : 1993 : Total Size= 8520 bytes File Size = 146 * *Baskets : 1 : Basket Size= 32000 bytes Compression= 55.08 * *............................................................................* *Br 3 :type : type/I * *Entries : 1993 : Total Size= 8510 bytes File Size = 145 * *Baskets : 1 : Basket Size= 32000 bytes Compression= 55.45 * *............................................................................* *Br 4 :r : r/F * *Entries : 1993 : Total Size= 8495 bytes File Size = 141 * *Baskets : 1 : Basket Size= 32000 bytes Compression= 57.00 * *............................................................................*The type encode which hypothesis the value of q corresponds to (observation has type zero). To see that, run:
%CODE{"cpp"}%
q->Draw("type")
gPad->SetLogy()
%ENDCODE%
![]() q needs to be multiplied by -2 in order to obtain the test-statistic -2*ln(L(A)/L(N)).
![]() ![]() Comments |