First, setup the environment with the Higgs Combination tools: https://twiki.cern.ch/twiki/bin/view/CMS/SWGuideHiggsAnalysisCombinedLimit
Currently working with 74X (check latest on HiggsCombine twiki).
(Taken from combine twiki on June 27)
export SCRAM_ARCH=slc6_amd64_gcc491
cmsrel CMSSW_7_4_7
cd CMSSW_7_4_7/src
cmsenv
git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit
cd HiggsAnalysis/CombinedLimit
git fetch origin
git checkout v6.3.1
scramv1 b clean
scramv1 b
(Needed for analytical reweighting. When trying to reproduce EPS results, check with Konstantin and Alexandra about which tag to use for the following repositories.)
cd ${CMSSW_BASE}/src/
git clone git@github.com:cms-hh/HHStatAnalysis.git
scramv1 b HHStatAnalysis/AnalyticalModels
cd ${CMSSW_BASE}/src/HHStatAnalysis
git clone git@github.com:cms-hh/Support.git
cd ${CMSSW_BASE}/src/HiggsAnalysis/
git clone git@github.com:ResonantHbbHgg/bbggLimits.git
cd ${CMSSW_BASE}/src/HiggsAnalysis/bbggLimits/
scramv1 b
In order to make limit trees from all samples use these script:
makeAllTrees.py -x nonres [--NRW]
Make non-res shape benchmark points trees (MVA based with 350 M(HH) threshold):
makeAllTrees.py -x nonres -f LT_OutDir \
--doCatMVA --MVAHMC0 0.970 --MVAHMC1 0.600 --MVALMC0 0.985 --MVALMC1 0.600 --massNR 350 --LMLJBTC 0.55 --LMSJBTC 0.55
You can also provide the locations of the flat trees if they are not the ones hard-coded in
the script, via -s, -d options.
In order to make the trees from the single Higgs samples, use --doSMHiggs option, and don't run over data (-d 0):
makeAllTrees.py -x nonres -f LT_OutDir -d 0 --doSMHiggs --genDiPhotonFilter \
--doCatMVA --MVAHMC0 0.970 --MVAHMC1 0.600 --MVALMC0 0.985 --MVALMC1 0.600 --massNR 350 --LMLJBTC 0.55 --LMSJBTC 0.55
Before we can proceed further with setting the limits, some hadding needs to be done:
for m in LowMass HighMass; do hadd -f LT_OutDir_${m}/LT_output_bbHToGG_M-125_13TeV_amcatnlo.root LT_OutDir_${m}/LT_output_bbHToGG_M-125_4FS_yb*.root; done
for m in LowMass HighMass; do hadd -f LT_OutDir_${m}/LT_output_GluGluToHHTo2B2G_AllNodes.root LT_OutDir_${m}/LT_output_GluGluToHHTo2B2G_node_*.root; done
The C++ Loop code to produce the Limit Trees is located at src/bbggLTMaker.cc. In order to run it over a single tree use the python script scripts/pyLimitTreeMaker.py, which exists in the $PATH after scram build. To run it just do:
pyLimitTreeMaker.py -f fileName.root -o outDir
where fileName.root is a an input Flat tree to be run over, and
outDir is where the output trees will be created. The
makeAllTrees.py mentioned in the beginning utilizes the
pyLimitTreeMaker.py and runs it over many input files.
More options for the pyLimitTreeMaker.py can be specified. To see all of them look
directly in the code.
In the non-resonant search, to get the limits at any parameter point of 5D space, we need
to reweigh the signal sample to that point. To do that we have a script,
scripts/MakeARWTree.py. Have a look at it to understand what it does. Then, to simplify
the production of those reweighted trees, we have another script which does everything on
batch, scripts/ArwTreesOnLSF.py. For example, to make trees for kl scan run:
python scripts/ArwTreesOnLSF.py -t KL
Once the limit trees are produced, we would like to make the 2D fits in (m(gg), m(bb)) plane, for each category, and then run the limits.
The main functions to do the fits are implemented in src/bbgg2DFitter.cc. The python
scripts are needed to handle many different situations (resonant, non-resonant,
re-weighting to various benchmark points, etc.). In order to run just one limit you need
scripts/pyLimits.py. Minimal options for the SM point are:
pyLimits.py -f conf_default.json -o outDir --nodes SM
Important: one has to specify the location of the input limit trees in
conf_default.json file. The above command must be run on lxplus, if the input root
files are located on EOS.
The pyLimits.py script would call runFullChain() method which is implemented in
python/LimitsUtil.py. So in fact, the LimitsUtil.py script is
the base code which interacts with the functions in bbgg2DFitter.cc.
Using the --nodes SM option tells it to use the Limit Tree produced from a single SM MC
sample.
Alternatively, one can do the limit on the re-weighted samples of the merged non-resonant
samples. (For the SM point this allows to increase the statistics of the signal
sample. Such re-weighting was used for EPS17 results of 2016 data.)
Run it like so:
pyLimits.py -f conf_default.json -o outDir --analyticalRW
In case of problems it's useful to increase verbosity level with -v 1(2,3) option. In
this case the logs should be found in your outDir/logs and in the master log,
outDir/mainLog_date-time.log
We have another script to facilitate running the limit for benchmarks, kl and kl-kt scans:
python scripts/runLimitsOnLSF.py -f conf_default.json -t [JHEP, KL, KLKT] [-o outDir]
The above command should give you the limits identical to the ones on SVN.
Good luck!
Make background fit plots for m(gg) and m(jj) in all categories:
source scripts/MakeSMHHFullBkgPlots.sh LIMSDIR
where LIMSDIR is a directory with the limits output. Similarly, for the signal shape (SM point), run:
python scripts/MakeSigPlotSimple.py LIMSDIR
In order make the non-resonant benchmark limit plot:
python scripts/MakeBenchmarksPlot.py LIMSDIR
To get the kl scan plot:
python scripts/MakeKLambdaScan.py LIMSDIR
For kl-kt scan plot, we first need to gather the results of all limits in a text file, and then run the plotting script:
python scripts/MakeKLKTScanTxtList.py LIMSDIR [-s]
python scripts/MakeKLKTplot.py -l LIMSDIR/KLKT_Scan_List.txt
Here, -s option can be used if only limits for kt>0 are produced. In this case the
plot is simply drawn symmetrically over (0,0) point in (kl,kt) coordinates.
Similarly to the non-resonant case, we first need to create the limit trees with categorization and then run the limits. Below are the commands needed to make the limit trees.
For low masses:
for s in Radion BulkGraviton;
do echo ${s};
for m in 250 260 270 280 300 320 340 350 400 450 500 550 600;
do echo ${s} ${m};
python scripts/makeAllTrees.py -x res -d /eos/cms/store/group/phys_higgs/resonant_HH/RunII/FlatTrees/2016/Mar82018_ForPubli_RafStyle/Data/Hadd/ -s /eos/cms/store/group/phys_higgs/resonant_HH/RunII/FlatTrees/2016/Mar82018_ForPubli_RafStyle/Signal/Hadd/ -f LT_ResOutDirOne_ --doCatMVA --MVAHMC0 0.960 --MVAHMC1 0.700 --MVALMC0 0.960 --MVALMC1 0.700 --massNR 500 --LMLJBTC 0.0 --LMSJBTC 0.0 --resMass ${m} --resType ${s};
done;
doneFor high masses (some numbers are different in the parameters wrt low masses):
for s in Radion BulkGraviton;
do echo ${s};
for m in 500 550 600 650 700 750 800 900;
do echo ${s} ${m};
python scripts/makeAllTrees.py -x res -d /eos/cms/store/group/phys_higgs/resonant_HH/RunII/FlatTrees/2016/Mar82018_ForPubli_RafStyle/Data/Hadd/ -s /eos/cms/store/group/phys_higgs/resonant_HH/RunII/FlatTrees/2016/Mar82018_ForPubli_RafStyle/Signal/Hadd/ -f LT_ResOutDirToo_ --doCatMVA --MVAHMC0 0.500 --MVAHMC1 0.000 --MVALMC0 0.500 --MVALMC1 0.000 --massNR 500 --LMLJBTC 0.0 --LMSJBTC 0.0 --resMass ${m} --resType ${s};
done;
doneTo produce the limits, modify the input paths in conf_resonant_LM.json and
conf_resonant_HM.json files and then run locally (it's rather quick compared to the
non-resonant limits):
pyLimits.py -f conf_resonant_LM.json --mass all-LM -j4 -o LIMS_RES_LM --overwrite
pyLimits.py -f conf_resonant_HM.json --mass all-HM -j4 -o LIMS_RES_HM --overwrite
(notice we run with separate config files for high mass points and low mass points, because the input paths for the limits trees must be different in those two cases).
Once the limits are produced, make a plot with this script:
python scripts/MakeResPlot.py LIMS_RES_LM --hmFolder=LIMS_RES_TEST_HM --log --isAsymptotic -n MyPlotName_Radi
python scripts/MakeResPlot.py LIMS_RES_LM --hmFolder=LIMS_RES_TEST_HM --log --isAsymptotic -n MyPlotName_Grav --isGrav