This aims to show the complete code workflow in detail.
In red are some comments on things that should be added or improved in the code, but I haven't had the time to do it.
In orange are some comments on things that I think(Matias) should be good to do, but are not necessary or important.
The main goal is to show how to go from the NanoAOD tuples to the mass fits to obtain the significance and asymptotic limits.
cmsrel CMSSW_14_0_14
cd CMSSW_14_0_14/src
cmsenv
git clone git@github.com:LPC-HH/HToMuMu.git
cd HmmAnalysis
make -j4
This is the first step. The code takes the NanoAOD tuples and saves the events that pass the selection in the analyser tuples. The selections are:
- Requieres two opposite charge muons with pt > 20GeV, abs(eta) < 2.4 and pass the medium identidication and loose isolation.
- One of the muons must have pt > 26 and match a HLT trigger muon.
For the Jets, the requirements are:
- pt > 25 GeV.
- abd(eta) < 4.7
- Pass loose isolation.
This code also checks if some of the jets are B-tagged based on the Particle net selection shown here: wiki
To execute the code, you can :
For simulations
./bin/HmmAnalyzer runList.txt out.root mc F 2016
For data
./bin/HmmAnalyzer runList.txt out.root data T 2016
To run over all the data sets, look at condor/README.md to run it in condor jobs
By default, these tuples are saved in:
/eos/uscms/store/user/lpchmumu/$(USER)/analyzer_HiggsMuMu/(Data/MC_background/MC_signal)
In order to obtain the weights, it is necessary to hadd the analyzer tuples to calculate the Sum of the generated weights.
Check condor/README.md to see how to hadd all the datasets. You can do that in the terminal using the hadd command, but it is recommended to use the script because it will check if some job failed before hadd the files.
By default, the hadded tuple is saved
This code takes the analyzer tuples and generates a new set of tuples, one for each dataset.
These tuples save just the variables that are being used now in the analysis, and also calculate variables that were not in the analyzer tuple as:
- weight: Event weight defined as (gen_weight x lumi x XS x pu_weight)/Sum_gen_weights.
- is_(ggH/VBF)_category: This variable is equal to 1 if the events pass the ggH/VBF selections.
Run the code for data/simulation :
./bin/CreateTuple analyzer_output.root output_directory era dataset T(data)/F(simulation) T(signal)/F(Bkg-data)
To run over all the data sets, look at condor/README.md to run it in condor jobs.
If you want to run it over a small dataset, look at scripts/README.md
By default, these tuples are saved in:
/eos/uscms/store/user/lpchmumu/$(USER)/analyzer_HiggsMuMu/tuples/
After you have the tuples, you need to have the sample. For example, you need to hadd all the 2022 Data samples in just one file or all the diBoson backgrounds in just one file. To do that, run
cd ./root_io/tuples/
bash hadd_tuples_EOS.sh
bash hadd_tuples.sh
If you want to add a variable to the tuples, you just need to go to ./lib/CreateTuple.h add the variable to the CreateTuple class, and to the CreateTuple::setBranchesAddressesInput() method.
For the training, you will need a VPN to Fermilab, do that first because it can take some time. See how to do it in ./python/xgboost/README.md
You need to have the tuples before continuing with this.
BDT does not like big files, so we need to generate a tuples that contain just the events from that production mode and the training variables.
To this use src/SkimTuples_VBF.cc src/SkimTuplesggH.cc and the scripts ./scripts/run_skim_ggH.py ./scripts/run_skim_VBF.py
You can change the datasets by changing the lists eras, background_datasets, and signal_datasets in the script.
Run
python3 ./scripts/run_skim_ggH.py
python3 ./scripts/run_skim_VBF.py
Once we have the skim tuples, you must hadd to have one background tuple and one signal sample, you can just do
cd ./root_io/skim/(ggH/VBF)
bash hadd_skim.sh
Once you have this, you can run the BDT training, look at ./python/xboost/README.md to see how to do it.
The key variables of this code are in ./python/xgboost/train.py are.
- DRAW_CORRELATIONS: if True generate a 2D map with the variables correations. (It takes time, so False for default).
- DO_STANDARDIZATION: if True standardize the variables (subtract the mean and divide by the standard deviation).
- USE_WEIGHT: if True add a weight the samples (gen_w * pu_weight * XS * lumi)/ SumGenWeights.
- USE_WEIGHT_MASS_RES: if true, add the mass resolution to the regular weight. Just one of the use_weight variable as true.
- APPEND_VARIABLES: if true, create new tuples with the BDT score saved on them.
- USE_BSCONSTRAIN: if true, use the bsConstrain variables for the training
The code will generate one file containing the model (python/xgboost/models/model_ggH_Combined_BFull_SNottH.xgb), one file containing the roc curve (python/xgboost/roc/ggH_Combined_BFull_SNottH_roc.txt), and will add the BDT score to the tuples. The tuples with the BDT score added are saved in root_io/tuples/BDT_score/(ggH/VBF)/BFull_SNottH/
Again, we need to hadd these tuples. To do this, do
bash scripts/hadd_bdt.sh
$\color{orange}\text{Here is a hardcoded path to VBF, the variable is called path. Change where it says VBF to ggH if you want that channel. It would be good to add an argument instead of changing the variable name.}$
This will copy the hadd_template.sh file in the directories where the tuples are, so the next time you can just do.
cd root_io/tuples/BDT_score/(ggH/VBF)/BFull_SNottH/
bash hadd.sh
Once we have this, we can find the categories. The categories are selected by scanning the cut x < BDT_score < 1 and finding which value maximizes S/sqrt(B) (Signal events divided by the square root of the background events). After this, we repeat the process but looking at x < BDT_score < cat1. This process is done 2 times for ggH and 3 times for VBF.
Run3
cd plot_macros
python3 plot_BDT_Categories.py (ggH/VBF) Combined
This will generate a plot showing the BDT output and the categories (plots/(ggH/VBF)_category/BDT_categories/BDT_output_Combined_BFull_SNottH.pdf), plots showing the BDT scan (/plots/ggH_category/BDT_categories/cuts), and also will print the categories, like this.
Categories: [0.0, 0.19184615384615372, 0.44615384615384546, 1.0]
Once you have this list, you need to manually copy these numbers in the dictionary call bdt_selections in the file ./plot_macros/plot_sim_vs_data.py and into the dictionary call bdt_cuts in the file ./scripts/split_tuples.py.
Inside the code, there is a variable called use_bsConstrain. If you train the BDT model using bsConstrain variable, set it to True, if you don't want to use bsConstrain variables set it to False.
Mass fits code based on this code Tutorial
I just changed a few things to make the code work for Hmumu(Br, xS, Lumi, etc), but the variables are still called mass_hgg or things like that.
The code requires a specific format for the tuples. To generate the necessary tuples run (be you that you complete the previous steps):
python3 script/split_tuples.py signal (ggH/VBF)
python3 script/split_tuples.py data (ggH/VBF)
Inside the code, there is a variable called use_bsConstrain. If you train the BDT model using bsConstrain variable, set it to True, if you don't want to use bsConstrain variables set it to False.
Once you have the tuples set up, the mass_fit_repo using the commands that are in the README
Once you have cloned the code (don't forget to run source setup.sh each time you close your lpc session), you can execute all the mass fits just by running the script.
- ./run_fits.sh : 4 paths, lines 11, 18, 24, 30.
- ./Signal/config_tutorial_combined.py : line 9
- ./Signal/config_tutorial_combined_VBF.py : line 9
- ./Background/config_tutorial.py : line 7
- ./Background/config_tutorial_VBF.py : line 7
bash run_fits.sh
This will generate all the necessary fits and will save them in flashggFinalFit/Signal/outdir_packaged, pflashggFinalFit/Background/outdir_tutorial, and flashggFinalFit/Plots/plots_tutorial.
Also, will create Datacards in this path: flashggFinalFit/Datacard
To understand what the scripts do, please read the whole Tutorial mentioned before. Also, here is a list of all the commands individually in case you want to execute them one by one to debug something
To calculate the limits, first do:
cd Combine
mkdir -p Models/signal
mkdir -p Models/background
cp ../Signal/outdir_packaged/CMS-HGG_sigfit_packaged*.root Models/signal/
cp ../Background/outdir_tutorial/CMS-HGG_multipdf*.root Models/background/
cp ../Datacard/Datacard_tutorial*.txt .
This will move the root files with the fits and the datacards to the right directory.
To go from datacarts.txt to datacards.root do:
python3 RunText2Workspace.py --mode mu_inclusive --batch local --ext _tutorial_(category)
where categorie can be ggH,,ggHcat1, gHcat2, gHcat3, VBF, VBFcat1, VBFcat2, VBFcat3, VBFcat4.
Then obtain the significance by running.
combine -M Significance Datacard_tutorial_(category)_mu_inclusive.root -m 125 -t -1 --expectSignal 1
And for the asymptotic limits:
combine -M AsymptoticLimits Datacard_tutorial_(category)_mu_inclusive.root -m 125 -t -1 --run blind
MC vs data plots are generated using ./plot_macros/plot_sim_vs_data.py. Inside of these files, there are two lists called variables and eras that determine which plots will be generated, so you can modify them if you don't want to generate all the plots. You can execute the code by running.
cd plot_macros/
python3 plot_sim_vs_data.py
If you want to plot the events for just one channel (ggH or VBF) there is a dictionary called bdt_cuts, you must change the ggH or VBF entry to [] and run
python3 plot_sim_vs_data.py (ggH/VBF)
If you want to plot the events for the different categories, there is a dictionary called bdt_cuts. You must change the ggH or VBF entry to [0,cat1,cat2...,1] and run
python3 plot_sim_vs_data.py (ggH/VBF)
cd plot_macros/
python3 plot_sig_vs_bkg.py (ggH/VBF)