Skip to content

Commit

Permalink
Merge pull request #2 from lucastorterotot/MSSM_CPV_dev
Browse files Browse the repository at this point in the history 
CPV model implementation
  • Loading branch information
@ArturAkh
ArturAkh authored Jan 8, 2021
2 parents 7174c5e + f6f7694 commit d97dc93
Show file tree
Hide file tree
Showing 5 changed files with 597 additions and 15 deletions.
48 changes: 46 additions & 2 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,8 @@ done;
## Datacard creation

```bash
morph_parallel.py --output output --analysis mssm --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/control_region_categories.txt --variable mt_tot_puppi --additional_arguments "--auto_rebin=1 --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root" --parallel 5
morph_parallel.py --output output --analysis mssm --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/mssm_classic_categories.txt --variable mt_tot_puppi --additional_arguments "--auto_rebin=1 --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root" --parallel 5
morph_parallel.py --output output --analysis mssm --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/control_region_categories.txt --variable mt_tot_puppi --additional_arguments="--auto_rebin=1" --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root --parallel 5
morph_parallel.py --output output --analysis mssm --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/mssm_classic_categories.txt --variable mt_tot_puppi --additional_arguments="--auto_rebin=1" --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root --parallel 5

for era in 2016 2017 2018;
do
Expand Down Expand Up @@ -164,6 +164,7 @@ is compared with the SM prediction.
* high mass no-btag categories contain only bbH, bbA and ggH, ggA
* btag categories & common control regions contain the full BSM signal model
* `mssm_vs_sm_h125`: same as above, but using for ggh the templates from the SM prediction, which are reweighted to the yield predicted by the MSSM scenario.
* `mssm_vs_sm_CPV` : here the considered Higgs bosons are `H1` (SM-like), `H2` and `H3` instead of the usual `h` `H` and `A`. _IN PROGRESS_

## Datacard creation for `mssm_vs_sm_classic_h125`

Expand Down Expand Up @@ -193,12 +194,25 @@ morph_parallel.py --output output --analysis mssm_vs_sm_h125 --eras 2016,2017,20
mkdir -p output_mssm_vs_sm_h125/combined/cmb/; rsync -av --progress output_mssm_vs_sm_h125/201?/htt_*/* output_mssm_vs_sm_h125/combined/cmb/
```

## Datacard creation for `mssm_vs_sm_CPV`

```bash
morph_parallel.py --output output --analysis mssm_vs_sm_CPV --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/control_region_categories.txt --variable mt_tot_puppi --parallel 5 --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root
morph_parallel.py --output output --analysis mssm_vs_sm_CPV --eras 2016,2017,2018 --category_list ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/mssm_classic_categories.txt --variable mt_tot_puppi --parallel 5 --sm_gg_fractions ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root

mkdir -p output_mssm_vs_sm_CPV/combined/cmb/; rsync -av --progress output_mssm_vs_sm_CPV/201?/htt_*/* output_mssm_vs_sm_CPV/combined/cmb/
```

## Workspace creation (exemplary for `mssm_vs_sm_classic_h125`)

```bash
ulimit -s unlimited
combineTool.py -M T2W -o ws_mh125.root -P CombineHarvester.MSSMvsSMRun2Legacy.MSSMvsSM:MSSMvsSM --PO filePrefix=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/ --PO modelFile=13,Run2017,mh125_13.root --PO MSSM-NLO-Workspace=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3_mssm_mode.root -i output_mssm_vs_sm_classic_h125/combined/cmb/ --PO minTemplateMass=110.0 --PO maxTemplateMass=3200.0
```
_NB_ : for the `mssm_vs_sm_CPV` analysis, use `mh1125_CPV_13.root` instead of `mh125_13.root` and `--PO MSSM-NLO-Workspace=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root` instead of `--PO MSSM-NLO-Workspace=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3_mssm_mode.root`. The model building file to use is `CombineHarvester.MSSMvsSMRun2Legacy.MSSMvsSM_CPV` with the model named `MSSMvsSM_CPV`. You may also change `ws_mh125.root` to `ws_mH1125_CPV.root`. The command to run would then be
```
combineTool.py -M T2W -o ws_mH1125_CPV.root -P CombineHarvester.MSSMvsSMRun2Legacy.MSSMvsSM_CPV:MSSMvsSM_CPV --PO filePrefix=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/ --PO modelFile=13,Run2017,mh1125_CPV_13.root --PO MSSM-NLO-Workspace=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/higgs_pt_v3.root -i output_mssm_vs_sm_CPV/combined/cmb/ --PO minTemplateMass=110.0 --PO maxTemplateMass=3200.0
```

## Model-dependent CLs 95% limits for `mssm_vs_sm_classic_h125`

Expand Down Expand Up @@ -289,6 +303,36 @@ combineTool.py -M AsymptoticGrid ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2
plotLimitGrid.py asymptotic_grid.root --scenario-label="M_{h}^{125} scenario (H,A#rightarrow#tau#tau)" --output mssm_mh125_mssm_vs_sm_heavy --title-right="137 fb^{-1} (13 TeV)" --cms-sub="Own Work" --contours="exp-2,exp-1,exp0,exp+1,exp+2,obs" --model_file=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/mh125_13_fixed.root --y-range 2.0,60.0 --x-title "m_{A} [GeV]"
```

## Model-dependent CLs 95% limits for `mssm_vs_sm_CPV`

**Computing limits:**

```bash
mkdir -p calculation_mh125_mssm_vs_sm_CPV
cd calculation_mh125_mssm_vs_sm_CPV

ulimit -s unlimited
combineTool.py -M AsymptoticGrid ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/mssm_asymptotic_grid_CPV.json -d ${CMSSW_BASE}/src/output_mssm_vs_sm_CPV/combined/cmb/ws_mH1125_CPV.root --job-mode 'condor' --task-name 'mssm_mh125_mssm_vs_sm_CPV_1' --dry-run -v1 --cminDefaultMinimizerStrategy 0 --X-rtd MINIMIZER_analytic --cminDefaultMinimizerTolerance 0.01 --redefineSignalPOI x --setParameters r=1 --freezeParameters r

# After adaption of each shell script and condor configuration matching mattern condor_mssm_mh125_mssm_vs_sm_CPV_1.{sh,sub}, submit to batch system:
condor_submit condor_mssm_mh125_mssm_vs_sm_CPV_1.sub
```

**Collecting limits:**

Basically the same command as above, but with a different task name (in case some of the jobs have failed)

```bash
ulimit -s unlimited
combineTool.py -M AsymptoticGrid ${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/input/mssm_asymptotic_grid_CPV.json -d ${CMSSW_BASE}/src/output_mssm_vs_sm_CPV/combined/cmb/ws_mH1125_CPV.root --job-mode 'condor' --task-name 'mssm_mh125_mssm_vs_sm_CPV_2' --dry-run -v1 --cminDefaultMinimizerStrategy 0 --X-rtd MINIMIZER_analytic --cminDefaultMinimizerTolerance 0.01 --redefineSignalPOI x --setParameters r=1 --freezeParameters r
```

**Plotting limits:**

```bash
plotLimitGrid.py asymptotic_grid.root --scenario-label="M_{H_{1}}^{125} CPV scenario (H_{1},H_{2},H_{3}#rightarrow#tau#tau)" --output mssm_mh125_mssm_vs_sm_CPV --title-right="137 fb^{-1} (13 TeV)" --cms-sub="Own Work" --contours="exp-2,exp-1,exp0,exp+1,exp+2,obs" --model_file=${CMSSW_BASE}/src/CombineHarvester/MSSMvsSMRun2Legacy/data/mh1125_CPV_13.root --y-range 1.0,20.0 --x-title "m_{H^{#pm}} [GeV]" --mass_histogram m_H1
```

## Adaptions needed for job submission via htcondor

The following line in file [CombineToolBase.py](https://github.com/KIT-CMS/CombineHarvester/blob/master/CombineTools/python/combine/CombineToolBase.py#L33) from the CombineHarvester package needs to be
Expand Down
70 changes: 59 additions & 11 deletions bin/MorphingMSSMvsSM.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -88,7 +88,7 @@ int main(int argc, char **argv) {

vector<string> mass_susy_ggH({}), mass_susy_qqH({}), parser_bkgs({}), parser_bkgs_em({}), parser_sm_signals({}), parser_main_sm_signals({});

string analysis = "sm"; // "sm", "mssm", "mssm_vs_sm_classic", "mssm_vs_sm_classic_h125", "mssm_vs_sm_heavy", "mssm_vs_sm", "mssm_vs_sm_h125"
string analysis = "sm"; // "sm", "mssm", "mssm_vs_sm_classic", "mssm_vs_sm_classic_h125", "mssm_vs_sm_heavy", "mssm_vs_sm", "mssm_vs_sm_h125", "mssm_vs_sm_CPV"
int era = 2016; // 2016, 2017 or 2018
po::variables_map vm;
po::options_description config("configuration");
Expand Down Expand Up @@ -182,6 +182,11 @@ int main(int argc, char **argv) {
mssm_ggH_signals = {"ggH_t", "ggH_b", "ggH_i","ggA_t", "ggA_b", "ggA_i"};
mssm_bbH_signals = {"bbH", "bbA"};
}
else if(analysis == "mssm_vs_sm_CPV")
{
mssm_ggH_signals = {"ggH1_t", "ggH1_b", "ggH1_i", "ggH2_t", "ggH2_b", "ggH2_i", "ggH3_t", "ggH3_b", "ggH3_i"};
mssm_bbH_signals = {"bbH1", "bbH2", "bbH3"};
}
mssm_signals = ch::JoinStr({mssm_ggH_signals, mssm_bbH_signals});

bkgs = {"EMB", "ZL", "TTL", "VVL", "jetFakes", "ggHWW125", "qqHWW125", "WHWW125", "ZHWW125"};
Expand Down Expand Up @@ -239,14 +244,21 @@ int main(int argc, char **argv) {
RooRealVar mh("mh", "mh", 125., 90., 4000.);
mA.setConstant(true);

// Define MSSM CPV model-dependent mass parameters mH3, mH2, mH1
RooRealVar mH3("mH3", "mH3", 125., 90., 4000.);
RooRealVar mH2("mH2", "mH2", 125., 90., 4000.);
RooRealVar mH1("mH1", "mH1", 125., 90., 4000.);
RooRealVar mHp("mHp", "mHp", 125., 90., 4000.);
mHp.setConstant(true);

// Define MSSM model-independent mass parameter MH
RooRealVar MH("MH", "MH", 125., 90., 4000.);
MH.setConstant(true);

// Define categories
map<string, Categories> cats;
// STXS stage 0 categories (optimized on ggH and VBF)
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm_classic_h125"){
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_CPV"){
cats["et"] = {
{ 1, "et_wjets_control"},
{ 32, "et_nobtag_tightmt"},
Expand Down Expand Up @@ -446,7 +458,7 @@ int main(int argc, char **argv) {
if(analysis == "sm"){
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, main_sm_signals, cats[chn], true);
}
else if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_classic_h125"){
else if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_CPV"){
if(analysis != "mssm_vs_sm" && analysis != "mssm_vs_sm_h125")
{
if(analysis != "mssm_vs_sm_classic_h125"){
Expand Down Expand Up @@ -508,9 +520,11 @@ int main(int argc, char **argv) {
}
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, ch::JoinStr({main_sm_signals, sm_signals}), cats[chn], true);
}
if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" ){
if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_CPV"){
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, ch::JoinStr({main_sm_signals, sm_signals}), cats[chn], true);
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, {"qqh"}, cats[chn], true);
if(analysis != "mssm_vs_sm_CPV"){
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, {"qqh"}, cats[chn], true);
}
if(analysis == "mssm_vs_sm_classic_h125")
{
cb.AddProcesses({""}, {"htt"}, {era_tag}, {chn}, {"ggh"}, cats[chn], true);
Expand Down Expand Up @@ -540,11 +554,25 @@ int main(int argc, char **argv) {
cb.cp().channel({chn}).process(main_sm_signals).ExtractShapes(
input_file_base, "$BIN/$PROCESS$MASS", "$BIN/$PROCESS$MASS_$SYSTEMATIC");
}
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125"){
if(analysis != "mssm_vs_sm_h125" and analysis != "mssm_vs_sm_classic_h125"){
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_CPV"){
if(analysis != "mssm_vs_sm_h125" and analysis != "mssm_vs_sm_classic_h125" and analysis != "mssm_vs_sm_CPV"){
cb.cp().channel({chn}).process(mssm_ggH_signals).ExtractShapes(
input_file_base, "$BIN/$PROCESS_$MASS", "$BIN/$PROCESS_$MASS_$SYSTEMATIC");
}
else if(analysis == "mssm_vs_sm_CPV")
{
// use the ggH_t,b,i shapes for all H1, H2, H3 and the bbH shape for H1, H2, H3 as a starting point
cb.cp().channel({chn}).process({"ggH1_i", "ggH2_i", "ggH3_i"}).ExtractShapes(
input_file_base, "$BIN/ggH_i_$MASS", "$BIN/ggH_i_$MASS_$SYSTEMATIC");
cb.cp().channel({chn}).process({"ggH1_t", "ggH2_t", "ggH3_t"}).ExtractShapes(
input_file_base, "$BIN/ggH_t_$MASS", "$BIN/ggH_t_$MASS_$SYSTEMATIC");
cb.cp().channel({chn}).process({"ggH1_b", "ggH2_b", "ggH3_b"}).ExtractShapes(
input_file_base, "$BIN/ggH_b_$MASS", "$BIN/ggH_b_$MASS_$SYSTEMATIC");

// Included in ggH1_x ?
// cb.cp().channel({chn}).process({"ggh"}).ExtractShapes(
// input_file_base, "$BIN/ggH125$MASS", "$BIN/ggH125$MASS_$SYSTEMATIC");
}
else
{
cb.cp().channel({chn}).process({"ggH_i", "ggH_t", "ggH_b", "ggA_i", "ggA_t", "ggA_b"}).ExtractShapes(
Expand All @@ -555,7 +583,7 @@ int main(int argc, char **argv) {
cb.cp().channel({chn}).process(mssm_bbH_signals).ExtractShapes(
input_file_base, "$BIN/bbH_$MASS", "$BIN/bbH_$MASS_$SYSTEMATIC");
}
if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125"){
if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_CPV"){
cb.cp().channel({chn}).process(ch::JoinStr({sm_signals,main_sm_signals})).ExtractShapes(
input_file_base, "$BIN/$PROCESS$MASS", "$BIN/$PROCESS$MASS_$SYSTEMATIC");
cb.cp().channel({chn}).process({"qqh"}).ExtractShapes(
Expand Down Expand Up @@ -951,7 +979,7 @@ int main(int argc, char **argv) {
}
}
// Desired Asimov model: BG + Higgs. Since H->tautau treated all as signal (together with mssm !!!), need to retrieve the SM H->tautau shapes & add it to the asimov dataset
else if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125"){
else if(analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_CPV"){
auto sm_signal_shape = cb.cp().bin({b}).process(ch::JoinStr({sm_signals, main_sm_signals})).GetShape();
std::cout << "[INFO] Integral of SM HTT signal shape in bin " << b << ": " << sm_signal_shape.Integral() << "\n";
bool no_background = (background_shape.GetNbinsX() == 1 && background_shape.Integral() == 0.0);
Expand Down Expand Up @@ -1051,17 +1079,37 @@ int main(int argc, char **argv) {
fractions_sm.Close();
}

if(analysis == "mssm_vs_sm_CPV")
{
mass_var = {
{"ggH1_t", &mH1}, {"ggH1_b", &mH1}, {"ggH1_i", &mH1},
{"ggH2_t", &mH2}, {"ggH2_b", &mH2}, {"ggH2_i", &mH2},
{"ggH3_t", &mH3}, {"ggH3_b", &mH3}, {"ggH3_i", &mH3},
{"bbH1", &mH1},
{"bbH2", &mH2},
{"bbH3", &mH3}
};

process_norm_map = {
{"ggH1_t", "prenorm"}, {"ggH1_b", "prenorm"}, {"ggH1_i", "prenorm"},
{"ggH2_t", "prenorm"}, {"ggH2_b", "prenorm"}, {"ggH2_i", "prenorm"},
{"ggH3_t", "prenorm"}, {"ggH3_b", "prenorm"}, {"ggH3_i", "prenorm"},
{"bbH1", "norm"},
{"bbH2", "norm"},
{"bbH3", "norm"}
};
}

dout("[INFO] Prepare demo.");
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125")
if(analysis == "mssm" || analysis == "mssm_vs_sm_classic" || analysis == "mssm_vs_sm_classic_h125" || analysis == "mssm_vs_sm_heavy" || analysis == "mssm_vs_sm" || analysis == "mssm_vs_sm_h125" || analysis == "mssm_vs_sm_CPV")
{
//TFile morphing_demo(("htt_mssm_morphing_" + category+ "_" + era_tag + "_" + analysis + "_demo.root").c_str(), "RECREATE");

// Perform morphing
std::cout << "[INFO] Performing template morphing for mssm ggh and bbh.\n";
auto morphFactory = ch::CMSHistFuncFactory();
morphFactory.SetHorizontalMorphingVariable(mass_var);
morphFactory.Run(cb, ws, process_norm_map);
morphFactory.Run(cb, ws, process_norm_map); // buggy with CPV at this stage

if(analysis == "mssm"){
// Adding 'norm' terms into workspace according to desired signals
Expand Down
13 changes: 13 additions & 0 deletions input/mssm_asymptotic_grid_CPV.json
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
{
"opts" : "--singlePoint 1.0",
"POIs" : ["mHp", "tanb"],
"grids" : [
["130:600|10", "1:20|1", ""],
["600:900|5", "1:7|1", ""],
["600:900|5", "7:13|0.5", ""],
["600:900|5", "13:20|1", ""],
["900:1000|25", "1:20|1", ""],
["1000:1500|50", "1:20|1", ""]
],
"hist_binning" : [137, 130, 1500, 20, 1, 20]
}
Loading

0 comments on commit d97dc93

Please sign in to comment.