diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0efeb534b3557358152e64ec86e5a78413468b35..b615f0aa8e86e5a6dcc51dc0417527d29dacddf3 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -4,15 +4,13 @@ stages:
before_script:
- export CI_HOME="$( pwd )"
- - git submodule update --init --recursive
- - export SCRAM_ARCH="slc6_amd64_gcc530"
- - export CMSSW_VERSION="CMSSW_8_0_26_patch2"
+ - export SCRAM_ARCH="slc7_amd64_gcc700"
+ - export CMSSW_VERSION="CMSSW_10_2_12"
- source /cvmfs/cms.cern.ch/cmsset_default.sh
- cd ..
- scramv1 project CMSSW $CMSSW_VERSION
- cd $CMSSW_VERSION/src
- eval `scramv1 runtime -sh`
- - scram setup $CI_HOME/data/gsl.xml
- cd $CI_HOME
compile:
@@ -21,13 +19,28 @@ compile:
- cvmfs
script:
- cd $CMSSW_BASE/src
- - mkdir -p TTH
+ - mkdir -p TTH/
- cp -r $CI_HOME TTH/
- - source TTH/CommonClassifier/setup/install_mem.sh
- - source TTH/CommonClassifier/setup/install_recoLikelihood.sh
+ - git clone https://gitlab.cern.ch/Zurich_ttH/MEIntegratorStandalone.git -b 10_2_X TTH/MEIntegratorStandalone
+ - git clone https://gitlab.cern.ch/kit-cn-cms-public/RecoLikelihoodReconstruction.git TTH/RecoLikelihoodReconstruction
+ - mkdir -p $CMSSW_BASE/lib/$SCRAM_ARCH/
+ - cp -R TTH/MEIntegratorStandalone/libs/* $CMSSW_BASE/lib/$SCRAM_ARCH/
+ - scram setup lhapdf
+ - scram setup TTH/MEIntegratorStandalone/deps/gsl.xml
- scram b
- cd $CMSSW_BASE
- find src -maxdepth 3 -type d | grep -e "^src/.*/.*/\(interface\|data\|python\)" | tar -czf $CI_HOME/cmssw.tgz lib biglib bin --exclude="*.pyc" --files-from -
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- cmssw.tgz
@@ -42,9 +55,20 @@ BDT:
- cd $CMSSW_BASE
- tar -xzf $CI_HOME/cmssw.tgz
- cd src
- - scram b python
+ - scram b
- cd $CI_HOME
- bdt_test | tee -i bdt.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- bdt.txt
@@ -61,6 +85,17 @@ DLBDT:
- scram b python
- cd $CI_HOME
- DLbdt_test | tee -i dlbdt.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- dlbdt.txt
@@ -77,6 +112,17 @@ MEM:
- scram b python
- cd $CI_HOME
- mem_test | tee -i mem.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- mem.txt
@@ -93,6 +139,17 @@ MEM_DL:
- scram b python
- cd $CI_HOME
- mem_test_dl | tee -i mem.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- mem_dl.txt
@@ -106,9 +163,20 @@ DNN_SL:
- cd $CMSSW_BASE
- tar -xzf $CI_HOME/cmssw.tgz
- cd src
- - scram b python
+ - scram b
- cd $CI_HOME
- dnn_SL_test | tee -i dnn_SL.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- dnn_SL.txt
@@ -122,9 +190,20 @@ DNN_DL:
- cd $CMSSW_BASE
- tar -xzf $CI_HOME/cmssw.tgz
- cd src
- - scram b python
+ - scram b
- cd $CI_HOME
- dnn_DL_test | tee -i dnn_DL.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- dnn_DL.txt
@@ -141,6 +220,17 @@ RecoLikelihood:
- scram b python
- cd $CI_HOME
- RecoLikelihoodVariables_test | tee -i RecoLikelihoodVariables.txt
+ only:
+ changes:
+ - bin/*.{c,cc,cpp,py,h,sh,xml}
+ - crab/*.{c,cc,cpp,py,h,sh,xml}
+ - data/*.{c,cc,cpp,py,h,sh,xml}
+ - interface/*.{c,cc,cpp,py,h,sh,xml}
+ - plugins/*.{c,cc,cpp,py,h,sh,xml}
+ - python/*.{c,cc,cpp,py,h,sh,xml}
+ - root/*.{c,cc,cpp,py,h,sh,xml}
+ - src/*.{c,cc,cpp,py,h,sh,xml}
+ - test/*.{c,cc,cpp,py,h,sh,xml}
artifacts:
paths:
- RecoLikelihoodVariables.txt
diff --git a/README.md b/README.md
index f684caa7c8893ae4c31542c74019ed07e3702e2d..9d82d1a8780964d0afa822530d9ea87746312552 100644
--- a/README.md
+++ b/README.md
@@ -1,24 +1,28 @@
Common classifiers for ttH analysis
===================================
-
+This branch is for analyses using CMSSW_10_2_X
Setup
-----
-~~~
+```bash
cd $CMSSW_BASE/src
cmsenv
-git clone ssh://git@gitlab.cern.ch:7999/ttH/CommonClassifier.git TTH/CommonClassifier
-source TTH/CommonClassifier/setup/install_mem.sh
-source TTH/CommonClassifier/setup/install_recoLikelihood.sh
+git clone ssh://git@gitlab.cern.ch:7999/ttH/CommonClassifier.git TTH/CommonClassifier -b 10_2_X
+git clone https://gitlab.cern.ch/Zurich_ttH/MEIntegratorStandalone.git -b 10_2_X TTH/MEIntegratorStandalone
+git clone https://gitlab.cern.ch/kit-cn-cms-public/RecoLikelihoodReconstruction.git TTH/RecoLikelihoodReconstruction
+mkdir -p $CMSSW_BASE/lib/$SCRAM_ARCH/
+cp -R TTH/MEIntegratorStandalone/libs/* $CMSSW_BASE/lib/$SCRAM_ARCH/
+scram setup lhapdf
+scram setup TTH/MEIntegratorStandalone/deps/gsl.xml
scram b
mem_test
mem_test_dl
bdt_test
dnn_test
-~~~
+```
Usage
-----
diff --git a/bin/mem_test.cc b/bin/mem_test.cc
index 2dc0ad752e67b5ea5edf517a5b3c36e06da2bd7a..bafe0195804e6725bde47a623b146df83bdf16b3 100644
--- a/bin/mem_test.cc
+++ b/bin/mem_test.cc
@@ -84,7 +84,7 @@ jets_variations.push_back(jet_variations);
//Assume that event category doesn't change
vector <bool> changes_jet_category;
-for (int i=0; abs(i)<abs(jet_variations[0].size())+1; i++){
+for (unsigned int i=0; i<jet_variations[0].size()+1; i++){
changes_jet_category.push_back(false);
}
diff --git a/bin/mem_test_dl.cc b/bin/mem_test_dl.cc
index dccc08522e59219cd7de45454d1bd5793f7e96b7..13eaf9d39710c170cc740cfc6b6196859ec3c87b 100644
--- a/bin/mem_test_dl.cc
+++ b/bin/mem_test_dl.cc
@@ -79,7 +79,7 @@ for (unsigned int ijet=0; ijet<jets_p4_nominal.size(); ijet++) {
//Assume for now that it doesn't modify the event category
vector <bool> changes_jet_category;
-for (int i=0; abs(i)<abs(jet_variations[0].size()) + 1; i++){
+for (unsigned int i=0; i < jet_variations[0].size() + 1; i++){
changes_jet_category.push_back(false);
}
diff --git a/data/btag_pdfs_2016_deepFlavour.root b/data/btag_pdfs_2016_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..3f9348d87565821a1a0da8643cc123d15a67ae1a
Binary files /dev/null and b/data/btag_pdfs_2016_deepFlavour.root differ
diff --git a/data/btag_pdfs_2017_deepFlavour.root b/data/btag_pdfs_2017_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..b357770bddd8a00406494b20d576be04fc659fbb
Binary files /dev/null and b/data/btag_pdfs_2017_deepFlavour.root differ
diff --git a/data/btag_pdfs_2018_deepFlavour.root b/data/btag_pdfs_2018_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..b8f011d7eb39da22f1307a0904de059642c7e48d
Binary files /dev/null and b/data/btag_pdfs_2018_deepFlavour.root differ
diff --git a/data/gsl.xml b/data/gsl.xml
deleted file mode 100644
index 40699148679294e230e8790e397b0c9ba9b482fa..0000000000000000000000000000000000000000
--- a/data/gsl.xml
+++ /dev/null
@@ -1,12 +0,0 @@
-<tool name="gsl" version="2.2.1">
- <info url="http://www.gnu.org/software/gsl/gsl.html"/>
- <lib name="gsl"/>
- <lib name="gslcblas"/>
- <client>
- <environment name="GSL_BASE" default="/cvmfs/cms.cern.ch/slc6_amd64_gcc530/external/gsl/2.2.1"/>
- <environment name="LIBDIR" default="$GSL_BASE/lib"/>
- <environment name="INCLUDE" default="$GSL_BASE/include"/>
- </client>
- <runtime name="ROOT_INCLUDE_PATH" value="$INCLUDE" type="path"/>
- <use name="root_cxxdefaults"/>
-</tool>
diff --git a/data/transfer_2016_deepFlavour.root b/data/transfer_2016_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..125eb5a1d682b46b72ed3178ce0d74a041cc6678
Binary files /dev/null and b/data/transfer_2016_deepFlavour.root differ
diff --git a/data/transfer_2017_deepFlavour.root b/data/transfer_2017_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..f9c31e9b503dd81d438dee5ee2156281103b0eb0
Binary files /dev/null and b/data/transfer_2017_deepFlavour.root differ
diff --git a/data/transfer_2018_deepFlavour.root b/data/transfer_2018_deepFlavour.root
new file mode 100644
index 0000000000000000000000000000000000000000..d87aac28b15db2af103dec5d14b90d02f633fd4f
Binary files /dev/null and b/data/transfer_2018_deepFlavour.root differ
diff --git a/interface/AngularVariables.h b/interface/AngularVariables.h
index a9a444bcb2a13eafbe6c77249553749831d7c255..d98d23efea0e219d1d3643320e1655899cc9bd73 100644
--- a/interface/AngularVariables.h
+++ b/interface/AngularVariables.h
@@ -23,7 +23,11 @@ class AngularVariables{
TLorentzVector Lepton;
TLorentzVector hadb;
TLorentzVector lepb;
+ TLorentzVector hadw;
+ TLorentzVector lepw;
TLorentzVector mET;
+
+ // mass templates
const double top_mass_had = 165;
const double top_mass_lep = 168;
const double W_mass = 80.0;
@@ -37,4 +41,4 @@ class AngularVariables{
};
-#endif
\ No newline at end of file
+#endif
diff --git a/interface/Chi2Reconstruction.h b/interface/Chi2Reconstruction.h
new file mode 100644
index 0000000000000000000000000000000000000000..8f02073608bf818fb3ca2f31b66f71764c8a78e2
--- /dev/null
+++ b/interface/Chi2Reconstruction.h
@@ -0,0 +1,87 @@
+#ifndef _Chi2Reconstruction_h
+#define _Chi2Reconstruction_h
+
+#include <vector>
+#include <map>
+#include "TVector.h"
+#include "TLorentzVector.h"
+
+//using namespace std;
+
+class Chi2Reconstruction{
+
+ public:
+ Chi2Reconstruction(std::string tag, std::string bosonName, double templateMass_Boson, double templateSigma_Boson);
+
+ virtual ~Chi2Reconstruction();
+
+ void InitReconstruction(const std::vector<TLorentzVector>& jets, const std::vector<double>& csvs, double CSVMwp);
+ void ReconstructNeutrino(const TLorentzVector& lepton, TLorentzVector met);
+
+ void ReconstructTTXSystem(const TLorentzVector& lepton, const std::vector<TLorentzVector>& jets, const std::vector<double>& csvs, TLorentzVector met, double CSVMwp);
+
+ double GetDeltaPhi(double phi1, double phi2);
+
+ // top system
+ std::map<std::string,double> GetVariableMap_TopSystemAngles();
+ std::map<std::string,double> GetVariableMap_TopSystemObjects();
+ std::map<std::string,double> FillVariableMap_TopSystemAngles();
+ std::map<std::string,double> FillVariableMap_TopSystemObjects();
+
+ // higgs system
+ std::map<std::string,double> GetVariableMap_BosonSystemAngles();
+ std::map<std::string,double> GetVariableMap_BosonSystemObjects();
+ std::map<std::string,double> FillVariableMap_BosonSystemAngles();
+ std::map<std::string,double> FillVariableMap_BosonSystemObjects();
+
+ private:
+ // tags
+ std::string tag;
+ std::string bosonName;
+
+ // variables for reconstruction
+ double btagCut;
+ double metPz[2];
+ double chi;
+
+ int nJets;
+
+ bool reconstructBoson;
+ int allowedMistags;
+
+ // best candidates
+ TLorentzVector best_Boson;
+ TLorentzVector best_B1;
+ TLorentzVector best_B2;
+
+ TLorentzVector best_topLep;
+ TLorentzVector best_topHad;
+ TLorentzVector best_bHad;
+ TLorentzVector best_bLep;
+ TLorentzVector best_wHad;
+ TLorentzVector best_wLep;
+ TLorentzVector best_lepton;
+ TLorentzVector best_mET;
+
+ // best chi2 values
+ double minChi = 1000000.;
+ double best_chi2_topHad = 10000.;
+ double best_chi2_topLep = 10000.;
+ double best_chi2_wHad = 10000.;
+ double best_chi2_Boson = 10000.;
+
+ // mass templates
+ const double templateMass_W = 83.7;
+ const double templateMass_topHad = 173.1;
+ const double templateMass_topLep = 168; // not updated
+
+ const double templateSigma_W = 11.9;
+ const double templateSigma_topHad = 18.5;
+ const double templateSigma_topLep = 26; // not updated
+
+ // template masses determined by init
+ double templateMass_Boson;
+ double templateSigma_Boson;
+};
+
+#endif
diff --git a/interface/HypothesisCombinatorics.h b/interface/HypothesisCombinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b3a904aefa84b97d938f3e67f1cfb0e1e094fe7
--- /dev/null
+++ b/interface/HypothesisCombinatorics.h
@@ -0,0 +1,61 @@
+#ifndef TTH_HYPOTHESISCOMBINATORICS_H
+#define TTH_HYPOTHESISCOMBINATORICS_H
+
+#include <vector>
+#include <map>
+#include "TLorentzVector.h"
+#include "TString.h"
+#include "TObjArray.h"
+#include "TObjString.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+#include "TTH/CommonClassifier/interface/PermutationManager.h"
+// #include "TTH/CommonClassifier/interface/BDTJAClassifier.h"
+
+// class to evaluate lepton plus jets BDT set
+class HypothesisCombinatorics{
+
+public:
+
+ HypothesisCombinatorics(const std::string& bdt_input_vars){
+ if( !bdt_input_vars.empty() ){
+ FillVariableNameList(bdt_input_vars, optional_varlist);
+ }
+ }
+
+ virtual ~HypothesisCombinatorics();
+
+ // setting the public function purely virtual forces the user to declare this function in a daughter class
+ virtual std::map<std::string, float> GetBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+
+ std::vector<std::string> GetVariableNames() const;
+ float* GetAdress(std::string variablelabel){return mvars->GetAdress(variablelabel);}
+ void SetBtagWP(const double& in_btagWP);
+
+ bool IsCentral(const TLorentzVector& vec, const double& central_eta_cut = 2.4){
+ return abs(vec.Eta())<central_eta_cut;
+ }
+protected:
+ // virtual std::map<std::string, float> FindBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+ // const std::vector<TLorentzVector> &selectedJetP4,
+ // const std::vector<double> &selectedJetCSV,
+ // const TLorentzVector &metP4){return std::map<std::string, float>(); }
+
+ void FillVariableNameList(const TString& variables, std::vector<std::string>& vec_to_fill);
+ // std::unique_ptr<BDTJAClassifier> bdtclassifier = nullptr;
+ std::vector<std::string> optional_varlist = std::vector<std::string>();
+ double btagWP;
+ std::unique_ptr<TMVA::Reader> reader_th;
+ std::unique_ptr<PermutationManager> permutator;
+ std::unique_ptr<MVAvarsBase> mvars;
+ std::string bdtoutput_name;
+ std::string bdt_name;
+
+ std::vector<std::string> BDTlabels;
+ unsigned int minJets;
+};
+
+#endif
\ No newline at end of file
diff --git a/interface/MVAvars.h b/interface/MVAvars.h
new file mode 100644
index 0000000000000000000000000000000000000000..44a0b78f39be9ace60693e2ee661fee52cd0e188
--- /dev/null
+++ b/interface/MVAvars.h
@@ -0,0 +1,31 @@
+#ifndef TTH_MVAVARS_H
+#define TTH_MVAVARS_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+#include "TTH/CommonClassifier/interface/MEMClassifier.h"
+#include "TTH/CommonClassifier/interface/AngularVariables.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+// class to evaluate lepton plus jets BDT set
+class MVAvars : public MVAvarsBase
+{
+
+ public:
+ MVAvars();
+ ~MVAvars();
+
+ void FillMVAvarMap(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+
+ private:
+ CommonBDTvars bdtvar;
+ MEMClassifier mem;
+};
+
+#endif
diff --git a/interface/MVAvarsBase.h b/interface/MVAvarsBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..9c3e0e69939e378c47e78d647b7ece419031a938
--- /dev/null
+++ b/interface/MVAvarsBase.h
@@ -0,0 +1,81 @@
+#ifndef TTH_MVAVARSBASE_H
+#define TTH_MVAVARSBASE_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+
+// class to evaluate lepton plus jets BDT set
+class MVAvarsBase
+{
+
+ public:
+ MVAvarsBase();
+ virtual ~MVAvarsBase();
+
+ // Call this method to return the MVa vars, provide all necessary inputs. Jet CSV should be sorted the same way as jet p4.
+ // We could also write a class to contain the jet CSV and p4 information
+ std::map<std::string, float> GetMVAvars(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+ std::map<std::string, float> GetMVAvars( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+
+ // return the variable names and their values for the last evaluated event
+ std::map<std::string, float>& GetVariables();
+ float* GetAdress(std::string variablelabel){return &variableMap[variablelabel];}
+ void SetWP(double WP);
+
+ double GetWP(){
+ return btagMcut;
+ }
+ virtual void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+
+ virtual void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ virtual std::map<std::string, TLorentzVector> GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ virtual bool SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx);
+
+ virtual void SetTopisleptonic(bool isit){cout << "ERROR: dont call `SetTopisleptonic()` for non MVAvarsJABDTthw-objects!" << endl;}
+
+ float GetGlobalDefault(){return globalDefault;}
+ std::vector<std::string> GetRecolabels(){return Recolabels;}
+
+ float GetHT(const std::vector<TLorentzVector> &selectedJetP4);
+ TLorentzVector GetLeptonicW(const TLorentzVector &leptonP4,
+ const TLorentzVector &metP4);
+
+ bool JetIsCentral(const TLorentzVector &jet){return TMath::Abs(jet.Eta()) < 2.4 && jet.Pt() > 30;}
+ bool JetIsSelected(const TLorentzVector &jet){return JetIsCentral(jet) || jet.Pt() > 40;}
+ bool JetIsTagged(const TLorentzVector &jet, const double& jetCSV){return jetCSV > btagMcut && JetIsCentral(jet);}
+
+ void ResetVariableMap();
+
+ protected:
+ std::map<std::string, float> variableMap;
+ std::vector<std::string> Recolabels;
+
+ double btagMcut = -99;
+ float globalDefault = -999.;
+};
+#endif
\ No newline at end of file
diff --git a/interface/MVAvarsJABDTthq.h b/interface/MVAvarsJABDTthq.h
new file mode 100644
index 0000000000000000000000000000000000000000..12b237e868f9169373c71ea2c99a501bef7b704e
--- /dev/null
+++ b/interface/MVAvarsJABDTthq.h
@@ -0,0 +1,40 @@
+#ifndef TTH_MVAVARSJABDTTHQ_H
+#define TTH_MVAVARSJABDTTHQ_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+// define unique indizes
+enum tHqIndexes {tHq_hdau1_idx, tHq_hdau2_idx, tHq_btop_idx, tHq_ljet_idx};
+
+// class to provide a variable container for tHq jet assignment hypothesis testing
+class MVAvarsJABDTthq : public MVAvarsBase
+{
+
+ public:
+ MVAvarsJABDTthq();
+ ~MVAvarsJABDTthq();
+
+ void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ std::map<std::string, TLorentzVector> GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ bool SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx);
+
+
+ private:
+};
+#endif
\ No newline at end of file
diff --git a/interface/MVAvarsJABDTthw.h b/interface/MVAvarsJABDTthw.h
new file mode 100644
index 0000000000000000000000000000000000000000..a97abadaf3be0e3143789ecfc11ed311d235f401
--- /dev/null
+++ b/interface/MVAvarsJABDTthw.h
@@ -0,0 +1,53 @@
+#ifndef TTH_MVAVARSJABDTTHW_H
+#define TTH_MVAVARSJABDTTHW_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+// define unique indizes
+enum tHWIndexes {tHW_btop_idx, tHW_whaddau1_idx, tHW_whaddau2_idx, tHW_hdau1_idx, tHW_hdau2_idx};
+
+// class to provide a variable container for thw jet assignment hypothesis testing
+class MVAvarsJABDTthw : public MVAvarsBase
+{
+
+ public:
+ MVAvarsJABDTthw();
+ ~MVAvarsJABDTthw();
+
+ void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ std::map<std::string, float> GetMVAvars( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+ {
+ FillMVAvarMap(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4, jets_idx);
+ return variableMap;
+ }
+
+
+ void SetTopisleptonic(bool isit){topisleptonic = isit;}
+
+ std::map<std::string, TLorentzVector> GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ bool SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx);
+
+ private:
+ bool topisleptonic;
+};
+#endif
\ No newline at end of file
diff --git a/interface/MVAvarsJABDTttbar.h b/interface/MVAvarsJABDTttbar.h
new file mode 100644
index 0000000000000000000000000000000000000000..792d23e22047bf442bcbe099875fa7662137b2ce
--- /dev/null
+++ b/interface/MVAvarsJABDTttbar.h
@@ -0,0 +1,39 @@
+#ifndef TTH_MVAVARSJABDTTTBAR_H
+#define TTH_MVAVARSJABDTTTBAR_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+// define unique indizes
+enum ttbarIndexes {ttbar_btoplep_idx, ttbar_whaddau1_idx, ttbar_whaddau2_idx, ttbar_btophad_idx};
+
+// class to provide a variable container for ttbar jet assignment hypothesis testing
+class MVAvarsJABDTttbar : public MVAvarsBase
+{
+
+ public:
+ MVAvarsJABDTttbar();
+ ~MVAvarsJABDTttbar();
+
+ void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ std::map<std::string, TLorentzVector> GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ bool SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx);
+
+ private:
+};
+#endif
\ No newline at end of file
diff --git a/interface/MVAvarsJABDTtth.h b/interface/MVAvarsJABDTtth.h
new file mode 100644
index 0000000000000000000000000000000000000000..24b6f3f7f07966383ea621ec6005db006f375606
--- /dev/null
+++ b/interface/MVAvarsJABDTtth.h
@@ -0,0 +1,39 @@
+#ifndef TTH_MVAVARSJABDTTTH_H
+#define TTH_MVAVARSJABDTTTH_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/CommonBDTvars.h"
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+// define unique indizes
+enum ttHIndexes {ttH_btoplep_idx, ttH_whaddau1_idx, ttH_whaddau2_idx, ttH_btophad_idx, ttH_hdau1_idx, ttH_hdau2_idx};
+
+// class to provide a variable container for ttH jet assignment hypothesis testing
+class MVAvarsJABDTtth : public MVAvarsBase
+{
+
+ public:
+ MVAvarsJABDTtth();
+ ~MVAvarsJABDTtth();
+
+ void FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ std::map<std::string, TLorentzVector> GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx);
+
+ bool SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx);
+
+ private:
+};
+#endif
\ No newline at end of file
diff --git a/interface/PermutationManager.h b/interface/PermutationManager.h
new file mode 100644
index 0000000000000000000000000000000000000000..02003fdaf336efae2208947d22df542752761beb
--- /dev/null
+++ b/interface/PermutationManager.h
@@ -0,0 +1,41 @@
+#ifndef THANALYSIS_PERMUTATIONMANAGER_H
+#define THANALYSIS_PERMUTATIONMANAGER_H
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include "TMath.h"
+
+
+
+using namespace std;
+
+class PermutationManager
+{
+ //
+ // construction / destruction
+ //
+public:
+ PermutationManager(unsigned int N_idx, unsigned int Jets_Max);
+ virtual ~PermutationManager();
+
+
+public:
+ unsigned int get_Npermutations(unsigned int Njets);
+ void get_permutation(std::vector<int>* permutation, unsigned int Njets, unsigned int i);
+
+ void show();
+ void show(unsigned int Njets);
+
+
+
+
+private:
+ const unsigned int Nidx;
+ const unsigned int JetsMax;
+ std::vector<std::vector<std::vector<int>>> permutations;
+
+ unsigned int getIndex(unsigned int Njets);
+};
+
+#endif
diff --git a/interface/ReconstructedVars.h b/interface/ReconstructedVars.h
new file mode 100644
index 0000000000000000000000000000000000000000..b28e86b225e9ab0eff78fb2f3732db5613e052a2
--- /dev/null
+++ b/interface/ReconstructedVars.h
@@ -0,0 +1,44 @@
+#ifndef TTH_ReconstructedVars_H
+#define TTH_ReconstructedVars_H
+#include <vector>
+#include <map>
+#include <math.h>
+#include "TLorentzVector.h"
+#include "TMVA/Reader.h"
+#include "TTH/CommonClassifier/interface/Chi2Reconstruction.h"
+
+// class to evaluate lepton plus jets BDT set
+class ReconstructedVars
+{
+
+ public:
+ ReconstructedVars(bool reconstruct_ttH, bool reconstruct_ttZ);
+ ~ReconstructedVars();
+
+ // Call this method to return the Reconstruced vars, provide all necessary inputs. Jet CSV should be sorted the same way as jet p4.
+ // We could also write a class to contain the jet CSV and p4 information
+ std::map<std::string, double> GetReconstructedVars(
+ const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+
+ // return the variable names with default initialization
+ std::map<std::string, double> GetVariables();
+ void SetWP(double WP);
+
+ private:
+ void ResetVariableMap();
+
+ std::string category;
+ std::map<std::string, double> variableMap;
+ double btagMcut = -99;
+
+ bool reconstruct_ttH;
+ bool reconstruct_ttZ;
+
+ Chi2Reconstruction chi2reco_ttH;
+ Chi2Reconstruction chi2reco_ttZ;
+};
+
+#endif
diff --git a/interface/thqHypothesisCombinatorics.h b/interface/thqHypothesisCombinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..f40089071a6e048981c30276e2fe112a929a218a
--- /dev/null
+++ b/interface/thqHypothesisCombinatorics.h
@@ -0,0 +1,20 @@
+#ifndef TTH_THQHYPOTHESISCOMBINATORICS_H
+#define TTH_THQHYPOTHESISCOMBINATORICS_H
+
+#include <vector>
+#include <map>
+#include "TLorentzVector.h"
+#include "TTH/CommonClassifier/interface/HypothesisCombinatorics.h"
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTthq.h"
+
+
+class thqHypothesisCombinatorics : public HypothesisCombinatorics
+{
+ public:
+ thqHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring);
+ ~thqHypothesisCombinatorics();
+
+ private:
+};
+
+#endif
diff --git a/interface/thwHypothesisCombinatorics.h b/interface/thwHypothesisCombinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..78eaf0ec756ca850641f08745832cb5dfabeda20
--- /dev/null
+++ b/interface/thwHypothesisCombinatorics.h
@@ -0,0 +1,25 @@
+#ifndef TTH_THWHYPOTHESISCOMBINATORICS_H
+#define TTH_THWHYPOTHESISCOMBINATORICS_H
+
+#include <vector>
+#include <map>
+#include "TLorentzVector.h"
+#include "TTH/CommonClassifier/interface/HypothesisCombinatorics.h"
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTthw.h"
+
+
+class thwHypothesisCombinatorics : public HypothesisCombinatorics
+{
+ public:
+ thwHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring);
+ ~thwHypothesisCombinatorics();
+
+ std::map<std::string, float> GetBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4);
+
+ private:
+};
+
+#endif
\ No newline at end of file
diff --git a/interface/ttbarHypothesisCombinatorics.h b/interface/ttbarHypothesisCombinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a13720de4f5a1238110dd1b938f18a467026c56
--- /dev/null
+++ b/interface/ttbarHypothesisCombinatorics.h
@@ -0,0 +1,20 @@
+#ifndef TTH_TTBARHYPOTHESISCOMBINATORICS_H
+#define TTH_TTBARHYPOTHESISCOMBINATORICS_H
+
+#include <vector>
+#include <map>
+#include "TLorentzVector.h"
+#include "TTH/CommonClassifier/interface/HypothesisCombinatorics.h"
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTttbar.h"
+
+
+class ttbarHypothesisCombinatorics : public HypothesisCombinatorics
+{
+ public:
+ ttbarHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring);
+ ~ttbarHypothesisCombinatorics();
+
+ private:
+};
+
+#endif
\ No newline at end of file
diff --git a/interface/tthHypothesisCombinatorics.h b/interface/tthHypothesisCombinatorics.h
new file mode 100644
index 0000000000000000000000000000000000000000..77e25c47c39bed69a7bcd9fab22035c69030ca03
--- /dev/null
+++ b/interface/tthHypothesisCombinatorics.h
@@ -0,0 +1,20 @@
+#ifndef TTH_TTHHYPOTHESISCOMBINATORICS_H
+#define TTH_TTHHYPOTHESISCOMBINATORICS_H
+
+#include <vector>
+#include <map>
+#include "TLorentzVector.h"
+#include "TTH/CommonClassifier/interface/HypothesisCombinatorics.h"
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTtth.h"
+
+
+class tthHypothesisCombinatorics : public HypothesisCombinatorics
+{
+ public:
+ tthHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring);
+ ~tthHypothesisCombinatorics();
+
+ private:
+};
+
+#endif
\ No newline at end of file
diff --git a/setup/install_mem.sh b/setup/install_mem.sh
deleted file mode 100644
index 13dafe7f864682abbff47c6473330e853bee5e16..0000000000000000000000000000000000000000
--- a/setup/install_mem.sh
+++ /dev/null
@@ -1,30 +0,0 @@
-action() {
- if [ -e "${CMSSW_BASE}" ]; then
- local origin="$( pwd )"
-
- # source this script to install the MEM package
- cd $CMSSW_BASE
-
- # make the MEM install dir
- mkdir -p src/TTH
- cd src/TTH
-
- # get the MEM code
- if [ -d MEIntegratorStandalone ]; then
- echo "MEIntegratorStandalone dir already exists, repository isn't cloned again."
- else
- git clone https://gitlab.cern.ch/jpata/MEIntegratorStandalone.git MEIntegratorStandalone --branch v0.4
- fi
-
- # copy the OpenLoops ME libraries
- mkdir -p ../../lib/$SCRAM_ARCH/
- cp -R MEIntegratorStandalone/libs/* ../../lib/$SCRAM_ARCH/
- eval `scramv1 runtime -sh`
- ls
-
- scram setup MEIntegratorStandalone/deps/gsl.xml
-
- cd $origin
- fi
-}
-action "$@"
diff --git a/src/AngularVariables.cpp b/src/AngularVariables.cpp
index e87320f83b29931ec5506c7a137c06faf29d87e1..c771579b6e337c75590639619950790bc6143272 100644
--- a/src/AngularVariables.cpp
+++ b/src/AngularVariables.cpp
@@ -30,6 +30,8 @@ void AngularVariables::ReconstructTopSystem(const TLorentzVector& lepton, const
Lepton.SetPxPyPzE(0.,0.,0.,0.);
hadb.SetPxPyPzE(0.,0.,0.,0.);
lepb.SetPxPyPzE(0.,0.,0.,0.);
+ hadw.SetPxPyPzE(0.,0.,0.,0.);
+ lepw.SetPxPyPzE(0.,0.,0.,0.);
mET.SetPxPyPzE(0.,0.,0.,0.);
double metPz[2]={0.,0.};
@@ -73,25 +75,26 @@ void AngularVariables::ReconstructTopSystem(const TLorentzVector& lepton, const
if (radical < 0.0)
{
- imaginary=true;
+ imaginary=true;
}
if(imaginary)
{
- radical=0;
- /*
- double value=.001;
- while(true)
- {
- met.SetPxPyPzE(pfmet_px,pfmet_py,0.0,sqrt(pow(pfmet_px,2)+pow(pfmet_py,2))); //neutrino mass 0, pt = sqrt(px^2+py^2)
- // energyLep = lepton.E();
- a = (W_mass*W_mass/(2.0*energyLep)) + (lepton.Px()*met.Px() + lepton.Py()*met.Py())/energyLep;
- radical = (2.0*lepton.Pz()*a/energyLep)*(2.0*lepton.Pz()*a/energyLep);
- radical = radical - 4.0*(1.0 - (lepton.Pz()/energyLep)*(lepton.Pz()/energyLep))*(met.Px()*met.Px() + met.Py()*met.Py()- a*a);
- if(radical>=0)
- break;
- pfmet_px-=pfmet_px*value;
- pfmet_py-=pfmet_py*value;
- }*/
+ radical=0;
+ /*
+ double value=.001;
+ while(true)
+ {
+ met.SetPxPyPzE(pfmet_px,pfmet_py,0.0,sqrt(pow(pfmet_px,2)+pow(pfmet_py,2))); //neutrino mass 0, pt = sqrt(px^2+py^2)
+ // energyLep = lepton.E();
+ a = (W_mass*W_mass/(2.0*energyLep)) + (lepton.Px()*met.Px() + lepton.Py()*met.Py())/energyLep;
+ radical = (2.0*lepton.Pz()*a/energyLep)*(2.0*lepton.Pz()*a/energyLep);
+ radical = radical - 4.0*(1.0 - (lepton.Pz()/energyLep)*(lepton.Pz()/energyLep))*(met.Px()*met.Px() + met.Py()*met.Py()- a*a);
+ if(radical>=0) break;
+
+ pfmet_px-=pfmet_px*value;
+ pfmet_py-=pfmet_py*value;
+ }
+ */
}
metPz[0] = (lepton.Pz()*a/energyLep) + 0.5*sqrt(radical);
@@ -112,49 +115,73 @@ void AngularVariables::ReconstructTopSystem(const TLorentzVector& lepton, const
/////////////////////////////////////////// Loop over all jets, both Pz, calcaulte chi-square ///////////////////////////////////////////
TLorentzVector metNew;
- for( int ipznu=0; ipznu<2; ipznu++ ){
+
+ // neurino solutions
+ for( int ipznu=0; ipznu<2; ipznu++ )
+ {
+
metNew.SetXYZM(met.Px(),met.Py(),metPz[ipznu],0.0); //neutrino has mass 0
//with b-tag info
std::vector<TLorentzVector> not_b_tagged,b_tagged;
+
//fill not_b_tagged and b_tagged
- for( int i=0;i<nJets;i++ ){
-
+ for( int i=0;i<nJets;i++ )
+ {
if( (csvs[i]>btagCut && i!=ind_second_lowest_csv && i!=ind_lowest_csv) ) b_tagged.push_back(jets[i]);
else not_b_tagged.push_back(jets[i]);
-
}
+
//first make possible t_lep's with b-tagged jets (includes making W_lep)
- for( int i=0; i<int(b_tagged.size()); i++ ){
- TLorentzVector W_lep=metNew+lepton; //used for histogram drawing only
- TLorentzVector top_lep=metNew+lepton+b_tagged.at(i);
- chi_top_lep=pow((top_lep.M()-top_mass_lep)/sigma_lepTop,2);
+ for( int i=0; i<int(b_tagged.size()); i++ )
+ {
+ // Leptonic W
+ TLorentzVector W_lep = metNew+lepton;
+ // Leptonic Top
+ TLorentzVector top_lep = metNew+lepton+b_tagged.at(i);
+
+ // Chi2 for leptonic top
+ chi_top_lep = pow((top_lep.M()-top_mass_lep)/sigma_lepTop,2);
+
//next make possible W_had's with not b-tagged jets
- for( int j=0; j<int(not_b_tagged.size()); j++ ){
- for( int k=0; k<int(not_b_tagged.size()); k++ ){
- if( j!=k ){
- TLorentzVector W_had=not_b_tagged.at(j)+not_b_tagged.at(k);
- chi_W_had=pow((W_had.M()-W_mass)/sigma_hadW,2);
+ for( int j=0; j<int(not_b_tagged.size()); j++ )
+ for( int k=0; k<int(not_b_tagged.size()); k++ )
+ {
+ if( j == k ) continue;
+
+ // Hadronic W
+ TLorentzVector W_had = not_b_tagged.at(j)+not_b_tagged.at(k);
+
+ // Chi2 for hadronic W
+ chi_W_had = pow((W_had.M()-W_mass)/sigma_hadW,2);
+
//now make possible top_had's (using the W_had + some b-tagged jet)
- for( int l=0; l<int(b_tagged.size()); l++ ){
- if( l!=i ){
- TLorentzVector top_had=W_had+b_tagged.at(l);
- chi_top_had=pow((top_had.M()-top_mass_had)/sigma_hadTop,2);
- chi=chi_top_lep+chi_W_had+chi_top_had;
- //accept the lowest chi
- if( chi<minChi ){
- minChi=chi;
- lepb = b_tagged.at(i);
- hadb = b_tagged.at(l);
- leptop = top_lep;
- hadtop = top_had;
- Lepton = lepton;
- mET = metNew;
- }
- }
+ for( int l=0; l<int(b_tagged.size()); l++ )
+ {
+ if( l==i ) continue;
+
+ // Hadronic Top
+ TLorentzVector top_had = W_had+b_tagged.at(l);
+
+ // Chi2 for hadronic top
+ chi_top_had = pow((top_had.M()-top_mass_had)/sigma_hadTop,2);
+
+ // combined chi2
+ chi=chi_top_lep+chi_W_had+chi_top_had;
+
+ //accept the lowest chi
+ if( chi < minChi ){
+ minChi = chi;
+ lepw = W_lep;
+ hadw = W_had;
+ lepb = b_tagged.at(i);
+ hadb = b_tagged.at(l);
+ leptop = top_lep;
+ hadtop = top_had;
+ Lepton = lepton;
+ mET = metNew;
}
}
}
- }
}
}
@@ -184,4 +211,4 @@ std::map<TString,double> AngularVariables::GetAngularVariables() {
return angular_variables;
-}
\ No newline at end of file
+}
diff --git a/src/Chi2Reconstruction.cpp b/src/Chi2Reconstruction.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..18e9cd6626df17630dfffe6b617c469348755fd4
--- /dev/null
+++ b/src/Chi2Reconstruction.cpp
@@ -0,0 +1,616 @@
+#include "../interface/Chi2Reconstruction.h"
+#include <iostream>
+using namespace std;
+
+Chi2Reconstruction::Chi2Reconstruction(std::string t, std::string n, double mBoson, double sBoson):
+ tag(t), bosonName(n), templateMass_Boson(mBoson), templateSigma_Boson(sBoson) {}
+
+Chi2Reconstruction::~Chi2Reconstruction() {}
+
+void Chi2Reconstruction::InitReconstruction(const std::vector<TLorentzVector>& jets, const std::vector<double>& csvs, double CSVMwp) {
+
+ // init best candidates empty
+ best_Boson.SetPxPyPzE(0.,0.,0.,0.);
+
+ best_B1.SetPxPyPzE(0.,0.,0.,0.);
+ best_B2.SetPxPyPzE(0.,0.,0.,0.);
+
+ best_topLep.SetPxPyPzE(0.,0.,0.,0.);
+ best_topHad.SetPxPyPzE(0.,0.,0.,0.);
+
+ best_bHad.SetPxPyPzE(0.,0.,0.,0.);
+ best_bLep.SetPxPyPzE(0.,0.,0.,0.);
+
+ best_wHad.SetPxPyPzE(0.,0.,0.,0.);
+ best_wLep.SetPxPyPzE(0.,0.,0.,0.);
+
+ best_lepton.SetPxPyPzE(0.,0.,0.,0.);
+ best_mET.SetPxPyPzE(0.,0.,0.,0.);
+
+ // z components of MET
+ metPz[0] = 0.;
+ metPz[1] = 0.;
+
+ // chi2 default value
+ chi = 1e6;
+ minChi = 1e5;
+
+ // get number of tags and btags
+ btagCut = CSVMwp;
+ nJets = int(jets.size());
+ int nBtags = 0;
+ for (int i = 0; i < nJets; i++)
+ if (csvs[i] > btagCut) nBtags++;
+ int nUntags = nJets-nBtags;
+
+ // only try to reconstruct Higgs/Z if more than six jets found
+ reconstructBoson = false;
+ if (nJets >= 6)
+ reconstructBoson = true;
+
+ // figure out number of allowed mistags
+ allowedMistags = 0;
+ // allow one mistag if only 4 jets
+ if(nJets==4) allowedMistags += 1;
+
+ // allow as many mistags as needed to account for hadronic W
+ if(nUntags==0) allowedMistags += 2;
+ if(nUntags==1) allowedMistags += 1;
+
+ // allow as many mistags as needed to account for 4 bs
+ if(reconstructBoson)
+ {
+ if(nBtags==2) allowedMistags += 2;
+ if(nBtags==3) allowedMistags += 1;
+ }
+
+ //std::cout << std::endl;
+ //std::cout << "jets|tags " << nJets << "|" << nBtags << std::endl;
+ //std::cout << "allowedMistags: " << allowedMistags << std::endl;
+
+}
+
+
+
+///////////////////////////// Neutrino /////////////////////////////////////////////////////////
+void Chi2Reconstruction::ReconstructNeutrino(const TLorentzVector& lepton, TLorentzVector met)
+{
+
+ double energyLep = lepton.E();
+ double a = (templateMass_W*templateMass_W/(2.0*energyLep)) + (lepton.Px()*met.Px() + lepton.Py()*met.Py())/energyLep;
+
+ double radical = (2.0*lepton.Pz()*a/energyLep)*(2.0*lepton.Pz()*a/energyLep);
+ radical = radical - 4.0*(1.0 - (lepton.Pz()/energyLep)*(lepton.Pz()/energyLep))*(met.Px()*met.Px() + met.Py()*met.Py()- a*a);
+ // set radical to zero if smaller than zero
+ if (radical < 0.0)
+ radical = 0;
+
+ metPz[0] = (lepton.Pz()*a/energyLep) + 0.5*sqrt(radical);
+ metPz[0] = metPz[0] / (1.0 - (lepton.Pz()/energyLep)*(lepton.Pz()/energyLep));
+ metPz[1] = (lepton.Pz()*a/energyLep) - 0.5*sqrt(radical);
+ metPz[1] = metPz[1] / (1.0 - (lepton.Pz()/energyLep)*(lepton.Pz()/energyLep));
+
+
+}
+
+
+void Chi2Reconstruction::ReconstructTTXSystem(const TLorentzVector& lepton, const std::vector<TLorentzVector>& jets, const std::vector<double>& csvs, TLorentzVector met, double CSVMwp)
+{
+
+ // initializers
+ //std::cout << "starting chi2 reco ========================================" << std::endl;
+ InitReconstruction(jets, csvs, CSVMwp);
+ ReconstructNeutrino(lepton, met);
+ TLorentzVector metNew;
+
+ int mistags = 0;
+ int mistagged_topLepB = 0;
+ int mistagged_topHadB = 0;
+ int mistagged_topHadQ1 = 0;
+ int mistagged_topHadQ2 = 0;
+ int mistagged_BosonB1 = 0;
+ int mistagged_BosonB2 = 0;
+
+ double chi2_topLep = 10000;
+ double chi2_topHad = 10000;
+ double chi2_wHad = 10000;
+ double chi2_Boson = 10000;
+
+ // neurino solutions
+ //std::cout << "looping over neutrino solutions" << std::endl;
+ for( int ipznu=0; ipznu<2; ipznu++ )
+ {
+ // get met vector for neutrino solution
+ metNew.SetXYZM(met.Px(),met.Py(),metPz[ipznu],0.0);
+
+ // first make possible leptonic tops with b-tagged jets (includes making wLep)
+ for (int i = 0; i < nJets; i++)
+ {
+ // leptonic W
+ TLorentzVector wLep = metNew+lepton;
+ // leptonic Top
+ TLorentzVector topLep = metNew+lepton+jets.at(i);
+
+ // chi2 for leptonic top
+ chi2_topLep = pow((topLep.M()-templateMass_topLep)/templateSigma_topLep, 2);
+ mistagged_topLepB = 0;
+ if (csvs.at(i) < btagCut) mistagged_topLepB = 1;
+
+ // next make possible wHads
+ for (int j = 0; j < nJets; j++)
+ {
+ if ( j == i ) continue;
+ for (int k = 0; k < nJets; k++)
+ {
+ if( k == i || k == j ) continue;
+
+ // hadronic W
+ TLorentzVector wHad = jets.at(j)+jets.at(k);
+
+ // chi2 for hadronic W
+ chi2_wHad = pow((wHad.M()-templateMass_W)/templateSigma_W, 2);
+ mistagged_topHadQ1 = 0;
+ mistagged_topHadQ2 = 0;
+ if (csvs.at(j) > btagCut) mistagged_topHadQ1 = 1;
+ if (csvs.at(k) > btagCut) mistagged_topHadQ2 = 1;
+
+ // next make possible hadronic top
+ for (int l = 0; l < nJets; l++)
+ {
+ if( l == i || l == j || l == k ) continue;
+
+ // hadronic top
+ TLorentzVector topHad = wHad+jets.at(l);
+
+ // chi2 for hadronic top
+ chi2_topHad = pow((topHad.M()-templateMass_topHad)/templateSigma_topHad, 2);
+ mistagged_topHadB = 0;
+ if (csvs.at(l) < btagCut) mistagged_topHadB = 1;
+
+ // count mistags so far
+ mistags = mistagged_topHadB+mistagged_topLepB+mistagged_topHadQ1+mistagged_topHadQ2;
+ if (mistags > allowedMistags) continue;
+
+ bool bosonReconstructed = false;
+ TLorentzVector boson;
+ TLorentzVector B1;
+ TLorentzVector B2;
+ boson.SetPxPyPzE(0.,0.,0.,0.);
+ B1.SetPxPyPzE(0.,0.,0.,0.);
+ B2.SetPxPyPzE(0.,0.,0.,0.);
+
+ if (reconstructBoson)
+ {
+ chi2_Boson = 1e6;
+ double current_chi2_Boson = 1e5;
+ double current_mistags;
+ TLorentzVector current_boson;
+
+ // next make possible higgs bosons
+ for (int m = 0; m < nJets; m++)
+ {
+ if( m == i || m == j || m == k || m == l ) continue;
+ for (int n = 0; n < nJets; n++)
+ {
+ if( n == i || n == j || n == k || n == l || n == m ) continue;
+
+ // Higgs
+ current_boson = jets.at(m)+jets.at(n);
+
+ // chi2 for higgs
+ current_chi2_Boson = pow((current_boson.M()-templateMass_Boson)/templateSigma_Boson, 2);
+ mistagged_BosonB1 = 0;
+ mistagged_BosonB2 = 0;
+ if (csvs.at(m) < btagCut) mistagged_BosonB1 = 1;
+ if (csvs.at(n) < btagCut) mistagged_BosonB2 = 1;
+
+ // count mistags
+ current_mistags = mistags+mistagged_BosonB1+mistagged_BosonB2;
+ if (current_mistags > allowedMistags) continue;
+
+ // save best solution
+ if(current_chi2_Boson < chi2_Boson)
+ {
+ bosonReconstructed = true;
+ mistags = current_mistags;
+ chi2_Boson = current_chi2_Boson;
+ boson = current_boson;
+ B1 = jets.at(m);
+ B2 = jets.at(n);
+ }
+
+ }
+ }
+ }
+
+ // continue if no higgs could be reconstructed
+ if (reconstructBoson && !bosonReconstructed) continue;
+
+ // combine chi2
+ chi = chi2_topLep + chi2_wHad + chi2_topHad;
+ if (reconstructBoson) chi+=chi2_Boson;
+
+ // accept lowest chi
+ if( chi < minChi )
+ {
+ minChi = chi;
+ best_chi2_topHad = chi2_topHad;
+ best_chi2_topLep = chi2_topLep;
+ best_chi2_wHad = chi2_wHad;
+ best_chi2_Boson = chi2_Boson;
+
+ best_wLep = wLep;
+ best_wHad = wHad;
+ best_bLep = jets.at(i);
+ best_bHad = jets.at(l);
+ best_topLep = topLep;
+ best_topHad = topHad;
+ best_lepton = lepton;
+ best_mET = metNew;
+ best_Boson = boson;
+ best_B1 = B1;
+ best_B2 = B2;
+ }
+ }
+ }
+ }
+ }
+
+
+ }
+ //if( nJets >= 6 )
+ //{
+ //std::cout << "================================" << std::endl;
+ //std::cout << "event printout for >=6jet events" << std::endl;
+ //std::cout << "allowed mistags " << allowedMistags << std::endl;
+ //std::cout << "minChi " << minChi << std::endl;
+ //std::cout << "bestChi_topHad " << best_chi2_topHad << std::endl;
+ //std::cout << "bestChi_topLep " << best_chi2_topLep << std::endl;
+ //std::cout << "bestChi_wHad " << best_chi2_wHad << std::endl;
+ //std::cout << "bestChi_Boson " << best_chi2_Boson << std::endl;
+ //std::cout << "hadtop pt " << best_topHad.Pt() << std::endl;
+ //std::cout << "higgs pt " << best_Boson.Pt() << std::endl;
+ //std::cout << "================================" << std::endl;
+ //}
+
+}
+
+
+
+///////////////////////// VARIABLE MAP INITIALIZATIONS ////////////////////////////////
+std::map<std::string, double> Chi2Reconstruction::GetVariableMap_BosonSystemAngles() {
+ std::map<std::string, double> varMap;
+
+ // Deta
+ varMap[tag+"Deta_"+bosonName+"_Lep"] = -999.;
+ varMap[tag+"Deta_"+bosonName+"_bLep"] = -999.;
+ varMap[tag+"Deta_"+bosonName+"_bHad"] = -999.;
+ varMap[tag+"Deta_"+bosonName+"_topHad"] = -999.;
+ varMap[tag+"Deta_"+bosonName+"_topLep"] = -999.;
+
+ // Dphi
+ varMap[tag+"Dphi_"+bosonName+"_Lep"] = -999.;
+ varMap[tag+"Dphi_"+bosonName+"_bLep"] = -999.;
+ varMap[tag+"Dphi_"+bosonName+"_bHad"] = -999.;
+ varMap[tag+"Dphi_"+bosonName+"_topHad"] = -999.;
+ varMap[tag+"Dphi_"+bosonName+"_topLep"] = -999.;
+
+ // cosdtheta
+ varMap[tag+"cosdTheta_"+bosonName+"_Lep"] = -999.;
+ varMap[tag+"cosdTheta_"+bosonName+"_bLep"] = -999.;
+ varMap[tag+"cosdTheta_"+bosonName+"_bHad"] = -999.;
+ varMap[tag+"cosdTheta_"+bosonName+"_topHad"] = -999.;
+ varMap[tag+"cosdTheta_"+bosonName+"_topLep"] = -999.;
+
+ return varMap;
+ }
+
+std::map<std::string, double> Chi2Reconstruction::GetVariableMap_TopSystemAngles() {
+ std::map<std::string, double> varMap;
+
+ // Deta
+ //varMap[tag+"Deta_Lep_bLep"] = -999.;
+ varMap[tag+"Deta_Lep_bHad"] = -999.;
+ varMap[tag+"Deta_Lep_topHad"] = -999.;
+ //varMap[tag+"Deta_Lep_wHad"] = -999.;
+ varMap[tag+"Deta_bLep_bHad"] = -999.;
+ varMap[tag+"Deta_bLep_topHad"] = -999.;
+ //varMap[tag+"Deta_bLep_wHad"] = -999.;
+ varMap[tag+"Deta_bLep_wLep"] = -999.;
+ varMap[tag+"Deta_topLep_bHad"] = -999.;
+ varMap[tag+"Deta_topLep_topHad"] = -999.;
+ //varMap[tag+"Deta_topLep_wHad"] = -999.;
+ //varMap[tag+"Deta_bHad_wLep"] = -999.;
+ varMap[tag+"Deta_bHad_wHad"] = -999.;
+ varMap[tag+"Deta_topHad_wLep"] = -999.;
+ //varMap[tag+"Deta_wHad_wLep"] = -999.;
+
+ // Dphi
+ //varMap[tag+"Dphi_Lep_bLep"] = -999.;
+ varMap[tag+"Dphi_Lep_bHad"] = -999.;
+ varMap[tag+"Dphi_Lep_topHad"] = -999.;
+ //varMap[tag+"Dphi_Lep_wHad"] = -999.;
+ varMap[tag+"Dphi_bLep_bHad"] = -999.;
+ varMap[tag+"Dphi_bLep_topHad"] = -999.;
+ //varMap[tag+"Dphi_bLep_wHad"] = -999.;
+ varMap[tag+"Dphi_bLep_wLep"] = -999.;
+ varMap[tag+"Dphi_topLep_bHad"] = -999.;
+ varMap[tag+"Dphi_topLep_topHad"] = -999.;
+ //varMap[tag+"Dphi_topLep_wHad"] = -999.;
+ //varMap[tag+"Dphi_bHad_wLep"] = -999.;
+ varMap[tag+"Dphi_bHad_wHad"] = -999.;
+ varMap[tag+"Dphi_topHad_wLep"] = -999.;
+ //varMap[tag+"Dphi_wHad_wLep"] = -999.;
+
+ // cos dtheta
+ //varMap[tag+"cosdTheta_Lep_bLep"] = -999.;
+ varMap[tag+"cosdTheta_Lep_bHad"] = -999.;
+ varMap[tag+"cosdTheta_Lep_topHad"] = -999.;
+ //varMap[tag+"cosdTheta_Lep_wHad"] = -999.;
+ varMap[tag+"cosdTheta_bLep_bHad"] = -999.;
+ varMap[tag+"cosdTheta_bLep_topHad"] = -999.;
+ //varMap[tag+"cosdTheta_bLep_wHad"] = -999.;
+ varMap[tag+"cosdTheta_bLep_wLep"] = -999.;
+ varMap[tag+"cosdTheta_topLep_bHad"] = -999.;
+ varMap[tag+"cosdTheta_topLep_topHad"] = -999.;
+ //varMap[tag+"cosdTheta_topLep_wHad"] = -999.;
+ //varMap[tag+"cosdTheta_bHad_wLep"] = -999.;
+ varMap[tag+"cosdTheta_bHad_wHad"] = -999.;
+ varMap[tag+"cosdTheta_topHad_wLep"] = -999.;
+ //varMap[tag+"cosdTheta_wHad_wLep"] = -999.;
+
+ return varMap;
+ }
+
+std::map<std::string, double> Chi2Reconstruction::GetVariableMap_BosonSystemObjects() {
+ std::map<std::string, double> varMap;
+
+ // higgs
+ varMap[tag+bosonName+"_M"] = -999.;
+ varMap[tag+bosonName+"_E"] = -999.;
+ varMap[tag+bosonName+"_Eta"] = -999.;
+ varMap[tag+bosonName+"_Phi"] = -999.;
+ varMap[tag+bosonName+"_Pt"] = -999.;
+
+ // b jets from higgs decay
+ varMap[tag+bosonName+"_BJet1_E"] = -999.;
+ varMap[tag+bosonName+"_BJet1_Eta"] = -999.;
+ varMap[tag+bosonName+"_BJet1_Phi"] = -999.;
+ varMap[tag+bosonName+"_BJet1_Pt"] = -999.;
+
+ varMap[tag+bosonName+"_BJet2_E"] = -999.;
+ varMap[tag+bosonName+"_BJet2_Eta"] = -999.;
+ varMap[tag+bosonName+"_BJet2_Phi"] = -999.;
+ varMap[tag+bosonName+"_BJet2_Pt"] = -999.;
+
+ // chi2
+ varMap[tag+"Chi2"+bosonName] = -999.;
+
+ return varMap;
+}
+
+std::map<std::string, double> Chi2Reconstruction::GetVariableMap_TopSystemObjects() {
+ std::map<std::string, double> varMap;
+
+ // hadronic top system
+ //varMap[tag+"TopHad_BJet_M"] = -999.;
+ varMap[tag+"TopHad_BJet_E"] = -999.;
+ varMap[tag+"TopHad_BJet_Eta"] = -999.;
+ varMap[tag+"TopHad_BJet_Phi"] = -999.;
+ varMap[tag+"TopHad_BJet_Pt"] = -999.;
+
+ //varMap[tag+"TopHad_W_M"] = -999.;
+ varMap[tag+"TopHad_W_E"] = -999.;
+ varMap[tag+"TopHad_W_Eta"] = -999.;
+ varMap[tag+"TopHad_W_Phi"] = -999.;
+ varMap[tag+"TopHad_W_Pt"] = -999.;
+
+ varMap[tag+"TopHad_M"] = -999.;
+ varMap[tag+"TopHad_E"] = -999.;
+ varMap[tag+"TopHad_Eta"] = -999.;
+ varMap[tag+"TopHad_Phi"] = -999.;
+ varMap[tag+"TopHad_Pt"] = -999.;
+
+
+ // leptonic top system
+ //varMap[tag+"TopLep_BJet_M"] = -999.;
+ varMap[tag+"TopLep_BJet_E"] = -999.;
+ varMap[tag+"TopLep_BJet_Eta"] = -999.;
+ varMap[tag+"TopLep_BJet_Phi"] = -999.;
+ varMap[tag+"TopLep_BJet_Pt"] = -999.;
+
+ //varMap[tag+"TopLep_W_M"] = -999.;
+ varMap[tag+"TopLep_W_E"] = -999.;
+ varMap[tag+"TopLep_W_Eta"] = -999.;
+ varMap[tag+"TopLep_W_Phi"] = -999.;
+ varMap[tag+"TopLep_W_Pt"] = -999.;
+
+ varMap[tag+"TopLep_M"] = -999.;
+ varMap[tag+"TopLep_E"] = -999.;
+ varMap[tag+"TopLep_Eta"] = -999.;
+ varMap[tag+"TopLep_Phi"] = -999.;
+ varMap[tag+"TopLep_Pt"] = -999.;
+
+ // chi2 values
+ varMap[tag+"Chi2Total"] = -999.;
+ varMap[tag+"Chi2TopHad"] = -999.;
+ varMap[tag+"Chi2TopLep"] = -999.;
+ varMap[tag+"Chi2WHad"] = -999.;
+
+ return varMap;
+ }
+
+
+
+double Chi2Reconstruction::GetDeltaPhi(double phi1, double phi2) {
+ double dphi = TMath::Abs( phi1-phi2 );
+ if(dphi >= TMath::Pi())
+ dphi = 2.*TMath::Pi() - dphi;
+
+ return dphi;
+ }
+
+///////////////////////// VARIABLE MAP FILLERS ////////////////////////////////
+std::map<std::string,double> Chi2Reconstruction::FillVariableMap_BosonSystemAngles() {
+ std::map<std::string, double> varMap;
+
+ // Deta
+ varMap[tag+"Deta_"+bosonName+"_Lep"] = TMath::Abs( best_Boson.Eta() - best_lepton.Eta() );
+ varMap[tag+"Deta_"+bosonName+"_bLep"] = TMath::Abs( best_Boson.Eta() - best_bLep.Eta() );
+ varMap[tag+"Deta_"+bosonName+"_bHad"] = TMath::Abs( best_Boson.Eta() - best_bHad.Eta() );
+ varMap[tag+"Deta_"+bosonName+"_topHad"] = TMath::Abs( best_Boson.Eta() - best_topHad.Eta() );
+ varMap[tag+"Deta_"+bosonName+"_topLep"] = TMath::Abs( best_Boson.Eta() - best_topLep.Eta() );
+
+ // Dphi
+ varMap[tag+"Dphi_"+bosonName+"_Lep"] = GetDeltaPhi( best_Boson.Phi(), best_lepton.Phi() );
+ varMap[tag+"Dphi_"+bosonName+"_bLep"] = GetDeltaPhi( best_Boson.Phi(), best_bLep.Phi() );
+ varMap[tag+"Dphi_"+bosonName+"_bHad"] = GetDeltaPhi( best_Boson.Phi(), best_bHad.Phi() );
+ varMap[tag+"Dphi_"+bosonName+"_topHad"] = GetDeltaPhi( best_Boson.Phi(), best_topHad.Phi() );
+ varMap[tag+"Dphi_"+bosonName+"_topLep"] = GetDeltaPhi( best_Boson.Phi(), best_topLep.Phi() );
+
+ // cosdtheta
+ varMap[tag+"cosdTheta_"+bosonName+"_Lep"] = TMath::Cos( (best_Boson.Vect()).Angle(best_lepton.Vect()) );
+ varMap[tag+"cosdTheta_"+bosonName+"_bLep"] = TMath::Cos( (best_Boson.Vect()).Angle(best_bLep.Vect()) );
+ varMap[tag+"cosdTheta_"+bosonName+"_bHad"] = TMath::Cos( (best_Boson.Vect()).Angle(best_bHad.Vect()) );
+ varMap[tag+"cosdTheta_"+bosonName+"_topHad"] = TMath::Cos( (best_Boson.Vect()).Angle(best_topHad.Vect()) );
+ varMap[tag+"cosdTheta_"+bosonName+"_topLep"] = TMath::Cos( (best_Boson.Vect()).Angle(best_topLep.Vect()) );
+
+
+ return varMap;
+}
+
+std::map<std::string,double> Chi2Reconstruction::FillVariableMap_TopSystemAngles() {
+ std::map<std::string, double> varMap;
+
+ // Deta
+ //varMap[tag+"Deta_Lep_bLep"] = TMath::Abs( best_lepton.Eta() - best_bLep.Eta() );
+ varMap[tag+"Deta_Lep_bHad"] = TMath::Abs( best_lepton.Eta() - best_bHad.Eta() );
+ varMap[tag+"Deta_Lep_topHad"] = TMath::Abs( best_lepton.Eta() - best_topHad.Eta() );
+ //varMap[tag+"Deta_Lep_wHad"] = TMath::Abs( best_lepton.Eta() - best_wHad.Eta() );
+ varMap[tag+"Deta_bLep_bHad"] = TMath::Abs( best_bLep.Eta() - best_bHad.Eta() );
+ varMap[tag+"Deta_bLep_topHad"] = TMath::Abs( best_bLep.Eta() - best_topHad.Eta() );
+ //varMap[tag+"Deta_bLep_wHad"] = TMath::Abs( best_bLep.Eta() - best_wHad.Eta() );
+ varMap[tag+"Deta_bLep_wLep"] = TMath::Abs( best_bLep.Eta() - best_wLep.Eta() );
+ varMap[tag+"Deta_topLep_bHad"] = TMath::Abs( best_topLep.Eta() - best_bHad.Eta() );
+ varMap[tag+"Deta_topLep_topHad"] = TMath::Abs( best_topLep.Eta() - best_topHad.Eta() );
+ //varMap[tag+"Deta_topLep_wHad"] = TMath::Abs( best_topLep.Eta() - best_wHad.Eta() );
+ //varMap[tag+"Deta_bHad_wLep"] = TMath::Abs( best_bHad.Eta() - best_wLep.Eta() );
+ varMap[tag+"Deta_bHad_wHad"] = TMath::Abs( best_bHad.Eta() - best_wHad.Eta() );
+ varMap[tag+"Deta_topHad_wLep"] = TMath::Abs( best_topHad.Eta() - best_wLep.Eta() );
+ //varMap[tag+"Deta_wHad_wLep"] = TMath::Abs( best_wHad.Eta() - best_wLep.Eta() );
+
+ // Dphi
+ //varMap[tag+"Dphi_Lep_bLep"] = GetDeltaPhi( best_lepton.Phi(), best_bLep.Phi() );
+ varMap[tag+"Dphi_Lep_bHad"] = GetDeltaPhi( best_lepton.Phi(), best_bHad.Phi() );
+ varMap[tag+"Dphi_Lep_topHad"] = GetDeltaPhi( best_lepton.Phi(), best_topHad.Phi() );
+ //varMap[tag+"Dphi_Lep_wHad"] = GetDeltaPhi( best_lepton.Phi(), best_wHad.Phi() );
+ varMap[tag+"Dphi_bLep_bHad"] = GetDeltaPhi( best_bLep.Phi(), best_bHad.Phi() );
+ varMap[tag+"Dphi_bLep_topHad"] = GetDeltaPhi( best_bLep.Phi(), best_topHad.Phi() );
+ //varMap[tag+"Dphi_bLep_wHad"] = GetDeltaPhi( best_bLep.Phi(), best_wHad.Phi() );
+ varMap[tag+"Dphi_bLep_wLep"] = GetDeltaPhi( best_bLep.Phi(), best_wLep.Phi() );
+ varMap[tag+"Dphi_topLep_bHad"] = GetDeltaPhi( best_topLep.Phi(), best_bHad.Phi() );
+ varMap[tag+"Dphi_topLep_topHad"] = GetDeltaPhi( best_topLep.Phi(), best_topHad.Phi() );
+ //varMap[tag+"Dphi_topLep_wHad"] = GetDeltaPhi( best_topLep.Phi(), best_wHad.Phi() );
+ //varMap[tag+"Dphi_bHad_wLep"] = GetDeltaPhi( best_bHad.Phi(), best_wLep.Phi() );
+ varMap[tag+"Dphi_bHad_wHad"] = GetDeltaPhi( best_bHad.Phi(), best_wHad.Phi() );
+ varMap[tag+"Dphi_topHad_wLep"] = GetDeltaPhi( best_topHad.Phi(), best_wLep.Phi() );
+ //varMap[tag+"Dphi_wHad_wLep"] = GetDeltaPhi( best_wHad.Phi(), best_wLep.Phi() );
+
+ // cos dtheta
+ //varMap[tag+"cosdTheta_Lep_bLep"] = TMath::Cos( (best_lepton.Vect()).Angle(best_bLep.Vect()) );
+ varMap[tag+"cosdTheta_Lep_bHad"] = TMath::Cos( (best_lepton.Vect()).Angle(best_bHad.Vect()) );
+ varMap[tag+"cosdTheta_Lep_topHad"] = TMath::Cos( (best_lepton.Vect()).Angle(best_topHad.Vect()) );
+ //varMap[tag+"cosdTheta_Lep_wHad"] = TMath::Cos( (best_lepton.Vect()).Angle(best_wHad.Vect()) );
+ varMap[tag+"cosdTheta_bLep_bHad"] = TMath::Cos( (best_bLep.Vect()).Angle(best_bHad.Vect()) );
+ varMap[tag+"cosdTheta_bLep_topHad"] = TMath::Cos( (best_bLep.Vect()).Angle(best_topHad.Vect()) );
+ //varMap[tag+"cosdTheta_bLep_wHad"] = TMath::Cos( (best_bLep.Vect()).Angle(best_wHad.Vect()) );
+ varMap[tag+"cosdTheta_bLep_wLep"] = TMath::Cos( (best_bLep.Vect()).Angle(best_wLep.Vect()) );
+ varMap[tag+"cosdTheta_topLep_bHad"] = TMath::Cos( (best_topLep.Vect()).Angle(best_bHad.Vect()) );
+ varMap[tag+"cosdTheta_topLep_topHad"] = TMath::Cos( (best_topLep.Vect()).Angle(best_topHad.Vect()) );
+ //varMap[tag+"cosdTheta_topLep_wHad"] = TMath::Cos( (best_topLep.Vect()).Angle(best_wHad.Vect()) );
+ //varMap[tag+"cosdTheta_bHad_wLep"] = TMath::Cos( (best_bHad.Vect()).Angle(best_wLep.Vect()) );
+ varMap[tag+"cosdTheta_bHad_wHad"] = TMath::Cos( (best_bHad.Vect()).Angle(best_wHad.Vect()) );
+ varMap[tag+"cosdTheta_topHad_wLep"] = TMath::Cos( (best_topHad.Vect()).Angle(best_wLep.Vect()) );
+ //varMap[tag+"cosdTheta_wHad_wLep"] = TMath::Cos( (best_wHad.Vect()).Angle(best_wLep.Vect()) );
+
+ return varMap;
+}
+
+std::map<std::string,double> Chi2Reconstruction::FillVariableMap_BosonSystemObjects() {
+ std::map<std::string, double> varMap;
+
+ // higgs
+ varMap[tag+bosonName+"_M"] = best_Boson.M();
+ varMap[tag+bosonName+"_E"] = best_Boson.E();
+ varMap[tag+bosonName+"_Eta"] = best_Boson.Eta();
+ varMap[tag+bosonName+"_Phi"] = best_Boson.Phi();
+ varMap[tag+bosonName+"_Pt"] = best_Boson.Pt();
+
+ // b jets from higgs decay
+ varMap[tag+bosonName+"_BJet1_E"] = best_B1.E();
+ varMap[tag+bosonName+"_BJet1_Eta"] = best_B1.Eta();
+ varMap[tag+bosonName+"_BJet1_Phi"] = best_B1.Phi();
+ varMap[tag+bosonName+"_BJet1_Pt"] = best_B1.Pt();
+
+ varMap[tag+bosonName+"_BJet2_E"] = best_B2.E();
+ varMap[tag+bosonName+"_BJet2_Eta"] = best_B2.Eta();
+ varMap[tag+bosonName+"_BJet2_Phi"] = best_B2.Phi();
+ varMap[tag+bosonName+"_BJet2_Pt"] = best_B2.Pt();
+
+ // chi2
+ varMap[tag+"Chi2"+bosonName] = best_chi2_Boson;
+
+ return varMap;
+}
+
+
+std::map<std::string,double> Chi2Reconstruction::FillVariableMap_TopSystemObjects() {
+ std::map<std::string, double> varMap;
+
+ // hadronic top system
+ //varMap[tag+"TopHad_BJet_M"] = best_bHad.M();
+ varMap[tag+"TopHad_BJet_E"] = best_bHad.E();
+ varMap[tag+"TopHad_BJet_Eta"] = best_bHad.Eta();
+ varMap[tag+"TopHad_BJet_Phi"] = best_bHad.Phi();
+ varMap[tag+"TopHad_BJet_Pt"] = best_bHad.Pt();
+
+ //varMap[tag+"TopHad_W_M"] = best_wHad.M();
+ varMap[tag+"TopHad_W_E"] = best_wHad.E();
+ varMap[tag+"TopHad_W_Eta"] = best_wHad.Eta();
+ varMap[tag+"TopHad_W_Phi"] = best_wHad.Phi();
+ varMap[tag+"TopHad_W_Pt"] = best_wHad.Pt();
+
+ varMap[tag+"TopHad_M"] = best_topHad.M();
+ varMap[tag+"TopHad_E"] = best_topHad.E();
+ varMap[tag+"TopHad_Eta"] = best_topHad.Eta();
+ varMap[tag+"TopHad_Phi"] = best_topHad.Phi();
+ varMap[tag+"TopHad_Pt"] = best_topHad.Pt();
+
+
+ // leptonic top system
+ //varMap[tag+"TopLep_BJet_M"] = best_bLep.M();
+ varMap[tag+"TopLep_BJet_E"] = best_bLep.E();
+ varMap[tag+"TopLep_BJet_Eta"] = best_bLep.Eta();
+ varMap[tag+"TopLep_BJet_Phi"] = best_bLep.Phi();
+ varMap[tag+"TopLep_BJet_Pt"] = best_bLep.Pt();
+
+ //varMap[tag+"TopLep_W_M"] = best_wLep.M();
+ varMap[tag+"TopLep_W_E"] = best_wLep.E();
+ varMap[tag+"TopLep_W_Eta"] = best_wLep.Eta();
+ varMap[tag+"TopLep_W_Phi"] = best_wLep.Phi();
+ varMap[tag+"TopLep_W_Pt"] = best_wLep.Pt();
+
+ varMap[tag+"TopLep_M"] = best_topLep.M();
+ varMap[tag+"TopLep_E"] = best_topLep.E();
+ varMap[tag+"TopLep_Eta"] = best_topLep.Eta();
+ varMap[tag+"TopLep_Phi"] = best_topLep.Phi();
+ varMap[tag+"TopLep_Pt"] = best_topLep.Pt();
+
+ // chi2 values
+ varMap[tag+"Chi2Total"] = minChi;
+ varMap[tag+"Chi2TopHad"] = best_chi2_topHad;
+ varMap[tag+"Chi2TopLep"] = best_chi2_topLep;
+ varMap[tag+"Chi2WHad"] = best_chi2_wHad;
+
+ return varMap;
+}
diff --git a/src/HypothesisCombinatorics.cpp b/src/HypothesisCombinatorics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..878b1206c456c19fd9b181b71468ec5ab951bcd8
--- /dev/null
+++ b/src/HypothesisCombinatorics.cpp
@@ -0,0 +1,125 @@
+#include "TTH/CommonClassifier/interface/HypothesisCombinatorics.h"
+
+
+// HypothesisCombinatorics::HypothesisCombinatorics(const std::string& optional_varstring){
+// if( !optional_varstring.empty() ){
+// FillVariableNameList(optional_varstring, optional_varlist);
+// }
+// }
+
+HypothesisCombinatorics::~HypothesisCombinatorics()
+{
+
+}
+
+// std::map<std::string, float> HypothesisCombinatorics::GetBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+// const std::vector<TLorentzVector> &selectedJetP4,
+// const std::vector<double> &selectedJetCSV,
+// const TLorentzVector &metP4) const
+// {
+// std::cout << "calling FindBestPermutation\n";
+// std::map<std::string, float> outputmap = FindBestPermutation(selectedLeptonP4,selectedJetP4,selectedJetCSV,metP4);
+// return outputmap;
+// }
+
+// virtual std::map<std::string, float> HypothesisCombinatorics::FindBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+// const std::vector<TLorentzVector> &selectedJetP4,
+// const std::vector<double> &selectedJetCSV,
+// const TLorentzVector &metP4) const
+// {
+// return std::map<std::string, float>();
+// }
+
+// std::vector<std::string> HypothesisCombinatorics::GetVariableNames() const{
+// if( bdtclassifier!= nullptr) return bdtclassifier->GetListOfVariables();
+// else return std::vector<std::string>();
+// }
+
+void HypothesisCombinatorics::FillVariableNameList(const TString& variables, std::vector<std::string>& vec_to_fill){
+ TString tosplit = "";
+ if(variables.EndsWith(".csv")){
+ ifstream in(variables.Data());
+ if(!in.is_open()) return;
+ std::cout << "DEBUG: reading from file " << variables.Data() << std::endl;
+ std::string line;
+ while(std::getline(in, line)){
+ if(!tosplit.EqualTo("")) tosplit += "," + line;
+ else tosplit = line;
+ }
+ }
+ else tosplit = variables;
+ TObjArray* tokens = tosplit.Tokenize(",");
+ for(int i = 0; i<tokens->GetEntries(); i++)
+ {
+ vec_to_fill.push_back( ((TObjString *)(tokens->At(i)))->String().Data() );
+ }
+}
+
+void HypothesisCombinatorics::SetBtagWP(const double& in_btagWP){
+ mvars->SetWP(in_btagWP);
+}
+
+std::vector<std::string> HypothesisCombinatorics::GetVariableNames() const{
+ std::vector<std::string> outvector;
+
+ for(auto& it : mvars->GetVariables()){
+ outvector.push_back(it.first);
+ }
+
+ outvector.push_back(bdtoutput_name);
+ return outvector;
+}
+
+std::map<std::string, float> HypothesisCombinatorics::GetBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ Float_t best_bdtout = -9999;
+ std::vector<int> permutation;
+ std::vector<int> best_reco_idx;
+ std::map<std::string, float> varmap;
+
+ //TODO jet cleaning
+
+ unsigned int NJets = selectedJetP4.size();
+
+ if(NJets < minJets)
+ {
+ // std::cout << "WARNING: found not enough jets! (found: " << NJets << ", required: " << minJets << ")\n";
+ mvars->ResetVariableMap();
+ }
+ else
+ {
+ // iterate permutations to find best
+ for(unsigned int permutation_idx=0; permutation_idx<permutator->get_Npermutations(NJets); permutation_idx++)
+ {
+ permutator->get_permutation(&permutation, NJets, permutation_idx);
+ if(mvars->SkipEvent(selectedJetP4, selectedJetCSV, permutation)) continue;
+
+ mvars->FillMVAvarMap(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4, permutation);
+
+ // eval BDT
+ Float_t recbdtout = reader_th->EvaluateMVA(bdt_name);
+
+ if(recbdtout > best_bdtout)
+ {
+ best_bdtout = recbdtout;
+ best_reco_idx = permutation;
+ }
+ }
+ if(best_reco_idx.size() == 0){
+// std::cout << "WARNING: found no list of best permutation indices! Resetting vars\n";
+ mvars->ResetVariableMap();
+ }
+ else
+ {
+ mvars->FillMVAvarMap(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4, best_reco_idx);
+ }
+ }
+
+ varmap = mvars->GetVariables();
+ varmap[bdtoutput_name.c_str()] = best_bdtout;
+
+ return varmap;
+}
\ No newline at end of file
diff --git a/src/MEMClassifier.cc b/src/MEMClassifier.cc
index 6c10cc8a193208b18b1d6b4303222dcf6ba24d2f..53279076c4f9e5dfb8a401812346e81c558ff0a2 100644
--- a/src/MEMClassifier.cc
+++ b/src/MEMClassifier.cc
@@ -354,7 +354,7 @@ std::vector<MEMResult> MEMClassifier::GetOutput(
vector<MEMResult> res;
- for (int i=0;abs(i)<abs(selectedJetVariation.size());i++) {
+ for (unsigned int i=0; i < selectedJetVariation.size();i++) {
MEMResult res2;
res.push_back(res2);
//Change hypothesis for the JES sources if needed
diff --git a/src/MVAvars.cpp b/src/MVAvars.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3314a3dc12ec9140497fc848ba010fe762704de
--- /dev/null
+++ b/src/MVAvars.cpp
@@ -0,0 +1,798 @@
+#include "TTH/CommonClassifier/interface/MVAvars.h"
+
+using namespace std;
+
+MVAvars::MVAvars()
+{
+ // ==================================================
+ //init all variables potentially used
+
+ // higgs specific variables
+ variableMap["Reco_best_higgs_mass"] = -999.;
+ variableMap["Reco_dEta_fn"] = -999.;
+
+
+ // CSV variables
+ variableMap["Evt_CSV_min"] = -999.;
+ variableMap["Evt_CSV_min_tagged"] = -999.;
+ variableMap["Evt_CSV_avg"] = -999.;
+ variableMap["Evt_CSV_avg_tagged"] = -999.;
+ variableMap["Evt_CSV_dev"] = -999.;
+ variableMap["Evt_CSV_dev_tagged"] = -999.;
+
+ // fox wolfram moments
+ variableMap["Evt_h0"] = -999.;
+ variableMap["Evt_h1"] = -999.;
+ variableMap["Evt_h2"] = -999.;
+ variableMap["Evt_h3"] = -999.;
+
+ // event masses
+ variableMap["Evt_M_Total"] = -999.;
+ variableMap["Evt_M3"] = -999.;
+ variableMap["Evt_M3_oneTagged"] = -999.;
+
+ // HT variables
+ variableMap["Evt_HT"] = -999.;
+ variableMap["Evt_HT_wo_MET"] = -999.;
+ variableMap["Evt_HT_jets"] = -999.;
+ variableMap["Evt_HT_tags"] = -999.;
+
+ // missing transveral variables
+ variableMap["Evt_MHT"] = -999.;
+ variableMap["Evt_MET"] = -999.;
+ variableMap["Evt_MTW"] = -999.;
+
+ // btag likelihood ratios
+ variableMap["Evt_blr"] = -999.;
+ variableMap["Evt_blr_transformed"] = -999.;
+
+ // Jet variables
+ variableMap["Evt_M_JetsAverage"] = -999.;
+ variableMap["Evt_Eta_JetsAverage"] = -999.;
+ variableMap["Evt_Pt_JetsAverage"] = -999.;
+ variableMap["Evt_E_JetsAverage"] = -999.;
+ variableMap["Evt_M_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_Eta_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_Pt_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_E_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_M_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_Eta_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_Pt_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_E_UntaggedJetsAverage"] = -999.;
+
+ // dijet variables (tagged jets)
+ variableMap["Evt_M2_closestTo125TaggedJets"] = -999.;
+ variableMap["Evt_M2_closestTo91TaggedJets"] = -999.;
+ variableMap["Evt_M2_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_Dr_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_Deta_TaggedJetsAverage"] = -999.;
+ variableMap["Evt_M2_minDrTaggedJets"] = -999.;
+ variableMap["Evt_Dr_minDrTaggedJets"] = -999.;
+ variableMap["Evt_Pt_minDrTaggedJets"] = -999.;
+ variableMap["Evt_Dr_maxDrTaggedJets"] = -999.;
+
+ // dijet variables (all jets)
+ variableMap["Evt_M2_JetsAverage"] = -999.;
+ variableMap["Evt_Dr_JetsAverage"] = -999.;
+ variableMap["Evt_Deta_JetsAverage"] = -999.;
+ variableMap["Evt_M2_minDrJets"] = -999.;
+ variableMap["Evt_Dr_minDrJets"] = -999.;
+ variableMap["Evt_Pt_minDrJets"] = -999.;
+ variableMap["Evt_Dr_maxDrJets"] = -999.;
+
+ // dijet variables (untagged jets)
+ variableMap["Evt_M2_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_Dr_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_Deta_UntaggedJetsAverage"] = -999.;
+ variableMap["Evt_M2_minDrUntaggedJets"] = -999.;
+ variableMap["Evt_Dr_minDrUntaggedJets"] = -999.;
+ variableMap["Evt_Pt_minDrUntaggedJets"] = -999.;
+ variableMap["Evt_Dr_maxDrUntaggedJets"] = -999.;
+
+ // max detas
+ variableMap["Evt_Deta_maxDetaJetJet"] = -999.;
+ variableMap["Evt_Deta_maxDetaTagTag"] = -999.;
+ variableMap["Evt_Deta_maxDetaJetTag"] = -999.;
+
+ // lepton + jet variables
+ variableMap["Evt_Dr_minDrLepTag"] = -999.;
+ variableMap["Evt_Dr_minDrLepJet"] = -999.;
+ variableMap["Evt_M_minDrLepTag"] = -999.;
+ variableMap["Evt_M_minDrLepJet"] = -999.;
+
+
+ // event shape variables
+ variableMap["Evt_aplanarity"] = -999.;
+ variableMap["Evt_aplanarity_jets"] = -999.;
+ variableMap["Evt_aplanarity_tags"] = -999.;
+ variableMap["Evt_sphericity"] = -999.;
+ variableMap["Evt_sphericity_jets"] = -999.;
+ variableMap["Evt_sphericity_tags"] = -999.;
+ variableMap["Evt_transverse_sphericity"] = -999.;
+ variableMap["Evt_transverse_sphericity_jets"] = -999.;
+ variableMap["Evt_transverse_sphericity_tags"] = -999.;
+ variableMap["Evt_JetPt_over_JetE"] = -999.;
+ variableMap["Evt_TaggedJetPt_over_TaggedJetE"] = -999.;
+
+ // reconstructed WLep variables
+ variableMap["Reco_WLep_E"] = -999.;
+ variableMap["Reco_WLep_Eta"] = -999.;
+ variableMap["Reco_WLep_Phi"] = -999.;
+ variableMap["Reco_WLep_Pt"] = -999.;
+ variableMap["Reco_WLep_Mass"] = -999.;
+
+
+}
+
+MVAvars::~MVAvars()
+{
+}
+
+void MVAvars::FillMVAvarMap(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ //check if CSVWP is set properly
+ if(btagMcut<0){
+ std::cerr << "Please set CSV Wp properly\n" << std::endl;
+ throw std::exception();
+ }
+ // Reset all map entries to globalDefault value so that noting is left over from the last event
+ ResetVariableMap();
+
+
+ // ==================================================
+ // construct object vectors etc
+
+ int njets = selectedJetP4.size();
+ int ntags = 0;
+ for (uint i = 0; i < selectedJetCSV.size(); i++)
+ {
+ if (selectedJetCSV[i] > btagMcut)
+ ntags++;
+ }
+
+ std::vector<double> selectedJetCSV_fixed;
+ std::vector<double> looseSelectedJetCSV;
+
+ for (uint i = 0; i < selectedJetCSV.size(); i++)
+ {
+ double tag = selectedJetCSV[i];
+ if (std::isnan(tag))
+ {
+ tag = -.1;
+ }
+ else if (tag < 0)
+ {
+ tag = -.1;
+ }
+ else if (tag > 1)
+ {
+ tag = 1.;
+ }
+ selectedJetCSV_fixed.push_back(tag);
+ // looseSelectedJetCSV.push_back(tag);
+ }
+
+ std::vector<TLorentzVector> looseSelectedJetP4;
+ std::vector<TLorentzVector> selectedTaggedJetP4;
+ std::vector<TLorentzVector> selectedUntaggedJetP4;
+ for (uint i = 0; i < selectedJetP4.size(); i++)
+ {
+ // looseSelectedJetP4.push_back(selectedJetP4[i]);
+ if (selectedJetCSV_fixed[i] > btagMcut)
+ {
+ selectedTaggedJetP4.push_back(selectedJetP4[i]);
+ }
+ else
+ {
+ selectedUntaggedJetP4.push_back(selectedJetP4[i]);
+ }
+ }
+
+ vector<vector<double>> jets_vvdouble;
+ for (auto jet = selectedJetP4.begin(); jet != selectedJetP4.end(); jet++)
+ {
+ vector<double> pxpypzE;
+ pxpypzE.push_back(jet->Px());
+ pxpypzE.push_back(jet->Py());
+ pxpypzE.push_back(jet->Pz());
+ pxpypzE.push_back(jet->E());
+ jets_vvdouble.push_back(pxpypzE);
+ }
+
+ vector<vector<double>> jets_vvdouble_tagged;
+ for (auto jet = selectedTaggedJetP4.begin(); jet != selectedTaggedJetP4.end(); jet++)
+ {
+ vector<double> pxpypzE;
+ pxpypzE.push_back(jet->Px());
+ pxpypzE.push_back(jet->Py());
+ pxpypzE.push_back(jet->Pz());
+ pxpypzE.push_back(jet->E());
+ jets_vvdouble_tagged.push_back(pxpypzE);
+ }
+
+ /// create jets_vvdouble_tagged
+ vector<double> sortedCSV = selectedJetCSV_fixed;
+ std::sort(sortedCSV.begin(), sortedCSV.end(), std::greater<double>());
+
+ // ==================================================
+ // calculate variables
+ // aplanarity and sphericity
+
+ double aplanarity, sphericity, transverse_sphericity;
+ bdtvar.getSp(selectedLeptonP4[0], metP4, selectedJetP4, aplanarity, sphericity, transverse_sphericity);
+
+ TLorentzVector dummy_Lepton;
+ TLorentzVector dummy_met;
+
+ double aplanarity_jets, sphericity_jets, transverse_sphericity_jets;
+ bdtvar.getSp(dummy_Lepton, dummy_met, selectedJetP4, aplanarity_jets, sphericity_jets, transverse_sphericity_jets);
+
+ double aplanarity_tags, sphericity_tags, transverse_sphericity_tags;
+ bdtvar.getSp(dummy_Lepton, dummy_met, selectedTaggedJetP4, aplanarity_tags, sphericity_tags, transverse_sphericity_tags);
+
+ // Fox Wolfram
+ double h0, h1, h2, h3, h4;
+ bdtvar.getFox(selectedJetP4, h0, h1, h2, h3, h4);
+
+ // best higgs mass 1
+ double minChi, dRbb;
+ TLorentzVector bjet1, bjet2;
+ double bestHiggsMass = -9.9;
+ // if (category == "6j4t")
+ if (njets >= 6 && ntags >= 4)
+ {
+ bestHiggsMass = bdtvar.getBestHiggsMass(selectedLeptonP4[0], metP4, selectedJetP4, selectedJetCSV_fixed, minChi, dRbb, bjet1, bjet2, looseSelectedJetP4, looseSelectedJetCSV, btagMcut);
+ }
+
+ // study top bb system
+ TLorentzVector dummy_metv;
+ double minChiStudy, chi2lepW, chi2leptop, chi2hadW, chi2hadtop, mass_lepW, mass_leptop, mass_hadW, mass_hadtop, dRbbStudy;
+ double testquant1, testquant2, testquant3, testquant4, testquant5, testquant6, testquant7;
+
+ TLorentzVector b1, b2;
+ bdtvar.study_tops_bb_syst(
+ metP4.Pt(), metP4.Phi(), dummy_metv, selectedLeptonP4[0], jets_vvdouble, selectedJetCSV_fixed,
+ minChiStudy, chi2lepW, chi2leptop, chi2hadW, chi2hadtop,
+ mass_lepW, mass_leptop, mass_hadW, mass_hadtop,
+ dRbbStudy,
+ testquant1, testquant2, testquant3, testquant4, testquant5, testquant6, testquant7,
+ b1, b2, btagMcut);
+
+ // fn
+ double dEta_fn = testquant6;
+ // ptE ratios
+ double pt_E_ratio = bdtvar.pt_E_ratio_jets(jets_vvdouble);
+ double pt_E_ratio_tag = bdtvar.pt_E_ratio_jets(jets_vvdouble_tagged);
+
+ // etamax
+ double jet_jet_etamax = bdtvar.get_jet_jet_etamax(jets_vvdouble);
+ double jet_tag_etamax = bdtvar.get_jet_tag_etamax(jets_vvdouble, selectedJetCSV_fixed, btagMcut);
+ double tag_tag_etamax = bdtvar.get_tag_tag_etamax(jets_vvdouble, selectedJetCSV_fixed, btagMcut);
+
+ // jet variables
+ double ht_jets = 0;
+ double ht_taggedjets = 0;
+ double Mlj = 0;
+ double dr_between_lep_and_closest_jet = 99;
+ double mht_px = 0;
+ double mht_py = 0;
+ TLorentzVector sum_pt_vec = selectedLeptonP4[0];
+ sum_pt_vec += metP4;
+ for (auto jetvec = selectedJetP4.begin(); jetvec != selectedJetP4.end(); ++jetvec)
+ {
+ ht_jets += jetvec->Pt();
+ mht_px += jetvec->Px();
+ mht_py += jetvec->Py();
+ sum_pt_vec += *jetvec;
+
+ double drLep = selectedLeptonP4[0].DeltaR(*jetvec);
+ if (drLep < dr_between_lep_and_closest_jet)
+ {
+ dr_between_lep_and_closest_jet = drLep;
+ Mlj = (selectedLeptonP4[0] + *jetvec).M();
+ }
+
+ }
+
+ mht_px += selectedLeptonP4[0].Px();
+ mht_py += selectedLeptonP4[0].Py();
+ double mass_of_everything = sum_pt_vec.M();
+ double ht_wo_met = ht_jets + selectedLeptonP4[0].Pt();
+ double ht = metP4.Pt() + ht_wo_met;
+ double MHT = sqrt(mht_px * mht_px + mht_py * mht_py);
+
+ double MTW = sqrt(2*(selectedLeptonP4[0].Pt()*metP4.Pt() - selectedLeptonP4[0].Px()*metP4.Px() - selectedLeptonP4[0].Py()*metP4.Py()));
+
+ // mass of lepton and closest bt-tagged jet
+ double Mlb = 0;
+ double dr_between_lep_and_closest_tagged_jet = 999.;
+ for (auto tagged_jet = selectedTaggedJetP4.begin(); tagged_jet != selectedTaggedJetP4.end(); tagged_jet++)
+ {
+ ht_taggedjets += tagged_jet->Pt();
+
+ double drLep = selectedLeptonP4[0].DeltaR(*tagged_jet);
+ if (drLep < dr_between_lep_and_closest_tagged_jet)
+ {
+ dr_between_lep_and_closest_tagged_jet = drLep;
+ Mlb = (selectedLeptonP4[0] + *tagged_jet).M();
+ }
+
+ }
+
+ // JET AVERAGES
+ // jet variables
+ double sumJetM = 0;
+ double sumJetPt = 0;
+ double sumJetE = 0;
+ double sumJetEta = 0;
+ int njets_all = 0;
+ for (auto itjet = selectedJetP4.begin(); itjet != selectedJetP4.end(); itjet++)
+ {
+ njets_all ++;
+ sumJetM += itjet->M();
+ sumJetPt += itjet->Pt();
+ sumJetE += itjet->E();
+ sumJetEta += itjet->Eta();
+ }
+ double avgJetM = -99.;
+ double avgJetPt = -99.;
+ double avgJetE = -99.;
+ double avgJetEta = -99.;
+ if (njets_all>0)
+ {
+ avgJetM = sumJetM/njets_all;
+ avgJetPt = sumJetPt/njets_all;
+ avgJetE = sumJetE/njets_all;
+ avgJetEta = sumJetEta/njets_all;
+ }
+
+ // tagged jet variables
+ double sumTaggedJetM = 0;
+ double sumTaggedJetPt = 0;
+ double sumTaggedJetE = 0;
+ double sumTaggedJetEta = 0;
+ int njets_tagged = 0;
+ for (auto itjet = selectedTaggedJetP4.begin(); itjet != selectedTaggedJetP4.end(); itjet++)
+ {
+ njets_tagged ++;
+ sumTaggedJetM += itjet->M();
+ sumTaggedJetPt += itjet->Pt();
+ sumTaggedJetE += itjet->E();
+ sumTaggedJetEta += itjet->Eta();
+ }
+ double avgTaggedJetM = -99.;
+ double avgTaggedJetPt = -99.;
+ double avgTaggedJetE = -99.;
+ double avgTaggedJetEta = -99.;
+ if (njets_tagged>0)
+ {
+ avgTaggedJetM = sumTaggedJetM/njets_tagged;
+ avgTaggedJetPt = sumTaggedJetPt/njets_tagged;
+ avgTaggedJetE = sumTaggedJetE/njets_tagged;
+ avgTaggedJetEta = sumTaggedJetEta/njets_tagged;
+ }
+
+ // untagged jet variables
+ double sumUntaggedJetM = 0;
+ double sumUntaggedJetPt = 0;
+ double sumUntaggedJetE = 0;
+ double sumUntaggedJetEta = 0;
+ int njets_untagged = 0;
+ for (auto itjet = selectedUntaggedJetP4.begin(); itjet != selectedUntaggedJetP4.end(); itjet++)
+ {
+ njets_untagged ++;
+ sumUntaggedJetM += itjet->M();
+ sumUntaggedJetPt += itjet->Pt();
+ sumUntaggedJetE += itjet->E();
+ sumUntaggedJetEta += itjet->Eta();
+ }
+ double avgUntaggedJetM = -99.;
+ double avgUntaggedJetPt = -99.;
+ double avgUntaggedJetE = -99.;
+ double avgUntaggedJetEta = -99.;
+ if (njets_untagged>0)
+ {
+ avgUntaggedJetM = sumUntaggedJetM/njets_untagged;
+ avgUntaggedJetPt = sumUntaggedJetPt/njets_untagged;
+ avgUntaggedJetE = sumUntaggedJetE/njets_untagged;
+ avgUntaggedJetEta = sumUntaggedJetEta/njets_untagged;
+ }
+
+
+ // DIJET SYSTEMS
+ // tagged jets
+ double closest_tagged_dijet_mass = -99;
+ double minDrTagged = 99;
+ double minPtTagged = -99;
+ double maxDrTagged = -99;
+ double sumDrTagged = 0;
+ double sumDetaTagged = 0;
+ double sumM2Tagged = 0;
+
+ int npairs_tagged = 0;
+ double tagged_dijet_mass_closest_to_125 = -999;
+ double tagged_dijet_mass_closest_to_91 = -999;
+ for (auto tagged_jet1 = selectedTaggedJetP4.begin(); tagged_jet1 != selectedTaggedJetP4.end(); tagged_jet1++)
+ {
+ for (auto tagged_jet2 = tagged_jet1 + 1; tagged_jet2 != selectedTaggedJetP4.end(); tagged_jet2++)
+ {
+ double dr = tagged_jet1->DeltaR(*tagged_jet2);
+ double m = (*tagged_jet1 + *tagged_jet2).M();
+ double pt = (*tagged_jet1 + *tagged_jet2).Pt();
+ sumDrTagged += dr;
+ sumM2Tagged += m;
+ double deta = fabs(tagged_jet1->Eta() - tagged_jet2->Eta());
+ sumDetaTagged += deta;
+ npairs_tagged++;
+ if (dr < minDrTagged)
+ {
+ minDrTagged = dr;
+ minPtTagged = pt;
+ closest_tagged_dijet_mass = m;
+ }
+ if (dr > maxDrTagged)
+ {
+ maxDrTagged = dr;
+ }
+ if (fabs(tagged_dijet_mass_closest_to_125 - 125) > fabs(m - 125))
+ {
+ tagged_dijet_mass_closest_to_125 = m;
+ }
+ if (fabs(tagged_dijet_mass_closest_to_91 - 91) > fabs(m - 91))
+ {
+ tagged_dijet_mass_closest_to_91 = m;
+ }
+ }
+ }
+
+ double avgDrTagged = -1;
+ double avgDetaTagged = 0;
+ double avgM2Tagged = 0;
+ if (npairs_tagged != 0)
+ {
+ avgDrTagged = sumDrTagged / (double)npairs_tagged;
+ avgDetaTagged = sumDetaTagged / (double)npairs_tagged;
+ avgM2Tagged = sumM2Tagged / (double)npairs_tagged;
+ }
+
+ // all jets
+ double closest_dijet_mass = -99;
+ double minDrJets = 99;
+ double minPtJets = -99;
+ double maxDrJets = -99;
+ double sumDrJets = 0;
+ double sumDetaJets = 0;
+ double sumM2Jets = 0;
+
+ int npairs_jets = 0;
+ for (auto jet1 = selectedJetP4.begin(); jet1 != selectedJetP4.end(); jet1++)
+ {
+ for (auto jet2 = jet1 + 1; jet2 != selectedJetP4.end(); jet2++)
+ {
+ double dr = jet1->DeltaR(*jet2);
+ double m = (*jet1 + *jet2).M();
+ double pt = (*jet1 + *jet2).Pt();
+ sumDrJets += dr;
+ sumM2Jets += m;
+ double deta = fabs(jet1->Eta() - jet2->Eta());
+ sumDetaJets += deta;
+ npairs_jets++;
+ if (dr < minDrJets)
+ {
+ minDrJets = dr;
+ minPtJets = pt;
+ closest_dijet_mass = m;
+ }
+ if (dr > maxDrJets)
+ {
+ maxDrJets = dr;
+ }
+ }
+ }
+
+ double avgDrJets = -1;
+ double avgDetaJets = 0;
+ double avgM2Jets = 0;
+ if (npairs_jets != 0)
+ {
+ avgDrJets = sumDrJets / (double)npairs_jets;
+ avgDetaJets = sumDetaJets / (double)npairs_jets;
+ avgM2Jets = sumM2Jets / (double)npairs_jets;
+ }
+
+ // untagged jets
+ double closest_untagged_dijet_mass = -99.;
+ double minDrUntagged = 99;
+ double minPtUntagged = -99;
+ double maxDrUntagged = -99;
+ double sumDrUntagged = 0;
+ double sumDetaUntagged = 0;
+ double sumM2Untagged = 0;
+
+ int npairs_untagged = 0;
+ for (auto jet1 = selectedUntaggedJetP4.begin(); jet1 != selectedUntaggedJetP4.end(); jet1++)
+ {
+ for (auto jet2 = jet1 + 1; jet2 != selectedUntaggedJetP4.end(); jet2++)
+ {
+ double dr = jet1->DeltaR(*jet2);
+ double m = (*jet1 + *jet2).M();
+ double pt = (*jet1 + *jet2).Pt();
+ sumDrUntagged += dr;
+ sumM2Untagged += m;
+ double deta = fabs(jet1->Eta() - jet2->Eta());
+ sumDetaUntagged += deta;
+ npairs_untagged++;
+ if (dr < minDrUntagged)
+ {
+ minDrUntagged = dr;
+ minPtUntagged = pt;
+ closest_untagged_dijet_mass = m;
+ }
+ if (dr > maxDrUntagged)
+ {
+ maxDrUntagged = dr;
+ }
+ }
+ }
+
+ double avgDrUntagged = -1;
+ double avgDetaUntagged = 0;
+ double avgM2Untagged = 0;
+ if (npairs_untagged != 0)
+ {
+ avgDrUntagged = sumDrUntagged / (double)npairs_untagged;
+ avgDetaUntagged = sumDetaUntagged / (double)npairs_untagged;
+ avgM2Untagged = sumM2Untagged / (double)npairs_untagged;
+ }
+
+
+ // M3
+ double m3 = -1.;
+ double maxpt_for_m3 = -1;
+ for (auto itJetVec1 = selectedJetP4.begin(); itJetVec1 != selectedJetP4.end(); ++itJetVec1)
+ {
+ for (auto itJetVec2 = itJetVec1 + 1; itJetVec2 != selectedJetP4.end(); ++itJetVec2)
+ {
+ for (auto itJetVec3 = itJetVec2 + 1; itJetVec3 != selectedJetP4.end(); ++itJetVec3)
+ {
+
+ TLorentzVector m3vec = *itJetVec1 + *itJetVec2 + *itJetVec3;
+
+ if (m3vec.Pt() > maxpt_for_m3)
+ {
+ maxpt_for_m3 = m3vec.Pt();
+ m3 = m3vec.M();
+ }
+ }
+ }
+ }
+
+ // M3 one tagged
+ double m3_onetagged = -1.;
+ double maxpt_for_m3_onetagged = -1.;
+ for (auto itTaggedJetVec1 = selectedTaggedJetP4.begin(); itTaggedJetVec1 != selectedTaggedJetP4.end(); ++itTaggedJetVec1)
+ {
+ for( auto itUntaggedJetVec1 = selectedUntaggedJetP4.begin(); itUntaggedJetVec1 != selectedUntaggedJetP4.end(); ++itUntaggedJetVec1)
+ {
+ for( auto itUntaggedJetVec2 = itUntaggedJetVec1+1; itUntaggedJetVec2 != selectedUntaggedJetP4.end(); ++itUntaggedJetVec2)
+ {
+ TLorentzVector m3vec_onetagged = *itTaggedJetVec1 + *itUntaggedJetVec1 + *itUntaggedJetVec2;
+ if (m3vec_onetagged.Pt() > maxpt_for_m3_onetagged)
+ {
+ maxpt_for_m3_onetagged = m3vec_onetagged.Pt();
+ m3_onetagged = m3vec_onetagged.M();
+ }
+ }
+ }
+ }
+
+ double detaJetsAverage = 0;
+ int nPairsJets = 0;
+ for (auto itJetVec1 = selectedJetP4.begin(); itJetVec1 != selectedJetP4.end(); ++itJetVec1)
+ {
+ for (auto itJetVec2 = itJetVec1 + 1; itJetVec2 != selectedJetP4.end(); ++itJetVec2)
+ {
+ detaJetsAverage += fabs(itJetVec1->Eta() - itJetVec2->Eta());
+ nPairsJets++;
+ }
+ }
+ if (nPairsJets > 0)
+ {
+ detaJetsAverage /= (double)nPairsJets;
+ }
+
+ // btag variables
+
+
+ // CSV averages
+ double averageCSV_tagged = 0;
+ double averageCSV_all = 0;
+ double lowest_btag_all = 99;
+ double lowest_btag_tagged = 99;
+ // int njets = selectedJetP4.size();
+ // int ntags = 0;
+ for (auto itCSV = selectedJetCSV_fixed.begin(); itCSV != selectedJetCSV_fixed.end(); ++itCSV)
+ {
+ averageCSV_all += fmax(*itCSV, 0);
+ lowest_btag_all = fmin(*itCSV, lowest_btag_all);
+ if (*itCSV < btagMcut)
+ continue;
+ lowest_btag_tagged = fmin(*itCSV, lowest_btag_tagged);
+ averageCSV_tagged += fmax(*itCSV, 0);
+ // ntags++;
+ }
+ if (ntags > 0)
+ averageCSV_tagged /= (double)ntags;
+ else
+ averageCSV_tagged = 0;
+ if (selectedJetCSV_fixed.size() > 0)
+ averageCSV_all /= selectedJetCSV_fixed.size();
+ else
+ averageCSV_all = 0;
+
+ if (lowest_btag_all > 90)
+ lowest_btag_all = -1;
+ if (lowest_btag_tagged > 90)
+ lowest_btag_tagged = -1;
+
+
+ // squared CSV deviations from average
+ double csvDev_all = 0;
+ double csvDev_tagged = 0;
+ for (auto itCSV = selectedJetCSV_fixed.begin(); itCSV != selectedJetCSV_fixed.end(); ++itCSV)
+ {
+ csvDev_all += pow(*itCSV - averageCSV_all, 2);
+ if (*itCSV < btagMcut)
+ continue;
+ csvDev_tagged += pow(*itCSV - averageCSV_tagged, 2);
+ }
+ if (ntags > 0)
+ csvDev_tagged /= (double)ntags;
+ else
+ csvDev_tagged = -1.;
+ if (selectedJetCSV_fixed.size() > 0)
+ csvDev_all /= selectedJetCSV_fixed.size();
+ else
+ csvDev_all = -1;
+
+
+
+ //calculate blr_transformed
+ double blr_transformed = -1.0;
+ std::vector<unsigned int> out_best_perm;
+ double out_P_4b = -1;
+ double out_P_2b = -1;
+ double eth_blr = -1;
+ if (selectedJetP4.size() > 3)
+ {
+ eth_blr = mem.GetBTagLikelihoodRatio(selectedJetP4,
+ selectedJetCSV_fixed,
+ out_best_perm,
+ out_P_4b,
+ out_P_2b);
+ }
+ blr_transformed = log(eth_blr / (1.0 - eth_blr));
+
+
+ // reconstruct Wlep
+ TLorentzVector Wlep = GetLeptonicW(selectedLeptonP4.at(0), metP4);
+
+ // ==================================================
+ // Fill variable map
+
+ // higgs specific variables
+ variableMap["Reco_best_higgs_mass"] = bestHiggsMass;
+ variableMap["Reco_dEta_fn"] = dEta_fn;
+
+ // CSV variables
+ variableMap["Evt_CSV_min"] = lowest_btag_all;
+ variableMap["Evt_CSV_min_tagged"] = lowest_btag_tagged;
+ variableMap["Evt_CSV_avg"] = averageCSV_all;
+ variableMap["Evt_CSV_avg_tagged"] = averageCSV_tagged;
+ variableMap["Evt_CSV_dev"] = csvDev_all;
+ variableMap["Evt_CSV_dev_tagged"] = csvDev_tagged;
+
+ // fox wolfram moments
+ variableMap["Evt_h0"] = h0;
+ variableMap["Evt_h1"] = h1;
+ variableMap["Evt_h2"] = h2;
+ variableMap["Evt_h3"] = h3;
+
+ // event masses
+ variableMap["Evt_M_Total"] = mass_of_everything;
+ variableMap["Evt_M3"] = m3;
+ variableMap["Evt_M3_oneTagged"] = m3_onetagged;
+
+ // HT variables
+ variableMap["Evt_HT"] = ht;
+ variableMap["Evt_HT_wo_MET"] = ht_wo_met;
+ variableMap["Evt_HT_jets"] = ht_jets;
+ variableMap["Evt_HT_tags"] = ht_taggedjets;
+
+ // missing transveral variables
+ variableMap["Evt_MHT"] = MHT;
+ variableMap["Evt_MET"] = metP4.Pt();
+ variableMap["Evt_MTW"] = MTW;
+
+ // btag likelihood ratios
+ variableMap["Evt_blr"] = eth_blr;
+ variableMap["Evt_blr_transformed"] = blr_transformed;
+
+ // Jet variables
+ variableMap["Evt_M_JetsAverage"] = avgJetM;
+ variableMap["Evt_Eta_JetsAverage"] = avgJetEta;
+ variableMap["Evt_Pt_JetsAverage"] = avgJetPt;
+ variableMap["Evt_E_JetsAverage"] = avgJetE;
+ variableMap["Evt_M_TaggedJetsAverage"] = avgTaggedJetM;
+ variableMap["Evt_Eta_TaggedJetsAverage"] = avgTaggedJetEta;
+ variableMap["Evt_Pt_TaggedJetsAverage"] = avgTaggedJetPt;
+ variableMap["Evt_E_TaggedJetsAverage"] = avgTaggedJetE;
+ variableMap["Evt_M_UntaggedJetsAverage"] = avgUntaggedJetM;
+ variableMap["Evt_Eta_UntaggedJetsAverage"] = avgUntaggedJetEta;
+ variableMap["Evt_Pt_UntaggedJetsAverage"] = avgUntaggedJetPt;
+ variableMap["Evt_E_UntaggedJetsAverage"] = avgUntaggedJetE;
+
+ // dijet variables (tagged jets)
+ variableMap["Evt_M2_closestTo125TaggedJets"] = tagged_dijet_mass_closest_to_125;
+ variableMap["Evt_M2_closestTo91TaggedJets"] = tagged_dijet_mass_closest_to_91;
+ variableMap["Evt_M2_TaggedJetsAverage"] = avgM2Tagged;
+ variableMap["Evt_Dr_TaggedJetsAverage"] = avgDrTagged;
+ variableMap["Evt_Deta_TaggedJetsAverage"] = avgDetaTagged;
+ variableMap["Evt_M2_minDrTaggedJets"] = closest_tagged_dijet_mass;
+ variableMap["Evt_Dr_minDrTaggedJets"] = minDrTagged;
+ variableMap["Evt_Pt_minDrTaggedJets"] = minPtTagged;
+ variableMap["Evt_Dr_maxDrTaggedJets"] = maxDrTagged;
+
+ // dijet variables (all jets)
+ variableMap["Evt_M2_JetsAverage"] = avgM2Jets;
+ variableMap["Evt_Dr_JetsAverage"] = avgDrJets;
+ variableMap["Evt_Deta_JetsAverage"] = avgDetaJets;
+ variableMap["Evt_M2_minDrJets"] = closest_dijet_mass;
+ variableMap["Evt_Dr_minDrJets"] = minDrJets;
+ variableMap["Evt_Pt_minDrJets"] = minPtJets;
+ variableMap["Evt_Dr_maxDrJets"] = maxDrJets;
+
+ // dijet variables (untagged jets)
+ variableMap["Evt_M2_UntaggedJetsAverage"] = avgM2Untagged;
+ variableMap["Evt_Dr_UntaggedJetsAverage"] = avgDrUntagged;
+ variableMap["Evt_Deta_UntaggedJetsAverage"] = avgDetaUntagged;
+ variableMap["Evt_M2_minDrUntaggedJets"] = closest_untagged_dijet_mass;
+ variableMap["Evt_Dr_minDrUntaggedJets"] = minDrUntagged;
+ variableMap["Evt_Pt_minDrUntaggedJets"] = minPtUntagged;
+ variableMap["Evt_Dr_maxDrUntaggedJets"] = maxDrUntagged;
+
+ // max detas
+ variableMap["Evt_Deta_maxDetaJetJet"] = jet_jet_etamax;
+ variableMap["Evt_Deta_maxDetaTagTag"] = tag_tag_etamax;
+ variableMap["Evt_Deta_maxDetaJetTag"] = jet_tag_etamax;
+
+ // lepton + jet variables
+ variableMap["Evt_Dr_minDrLepTag"] = dr_between_lep_and_closest_tagged_jet;
+ variableMap["Evt_Dr_minDrLepJet"] = dr_between_lep_and_closest_jet;
+ variableMap["Evt_M_minDrLepTag"] = Mlb;
+ variableMap["Evt_M_minDrLepJet"] = Mlj;
+
+ // event shape variables
+ variableMap["Evt_aplanarity"] = aplanarity;
+ variableMap["Evt_aplanarity_jets"] = aplanarity_jets;
+ variableMap["Evt_aplanarity_tags"] = aplanarity_tags;
+ variableMap["Evt_sphericity"] = sphericity;
+ variableMap["Evt_sphericity_jets"] = sphericity_jets;
+ variableMap["Evt_sphericity_tags"] = sphericity_tags;
+ variableMap["Evt_transverse_sphericity"] = transverse_sphericity;
+ variableMap["Evt_transverse_sphericity_jets"] = transverse_sphericity_jets;
+ variableMap["Evt_transverse_sphericity_tags"] = transverse_sphericity_tags;
+ variableMap["Evt_JetPt_over_JetE"] = pt_E_ratio;
+ variableMap["Evt_TaggedJetPt_over_TaggedJetE"] = pt_E_ratio_tag;
+
+ // reconstructed WLep variables
+ variableMap["Reco_WLep_E"] = Wlep.E();
+ variableMap["Reco_WLep_Eta"] = Wlep.Eta();
+ variableMap["Reco_WLep_Phi"] = Wlep.Phi();
+ variableMap["Reco_WLep_Pt"] = Wlep.Pt();
+ variableMap["Reco_WLep_Mass"] = Wlep.M();
+
+}
\ No newline at end of file
diff --git a/src/MVAvarsBase.cpp b/src/MVAvarsBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..91b69ec45acaee5c83e453f69f4290d513cec7c6
--- /dev/null
+++ b/src/MVAvarsBase.cpp
@@ -0,0 +1,143 @@
+#include "TTH/CommonClassifier/interface/MVAvarsBase.h"
+
+using namespace std;
+
+MVAvarsBase::MVAvarsBase():
+btagMcut(-99),
+globalDefault(-999.)
+{
+}
+
+MVAvarsBase::~MVAvarsBase()
+{
+}
+
+void MVAvarsBase::ResetVariableMap()
+{
+ for (auto it = variableMap.begin(); it != variableMap.end(); it++)
+ {
+ it->second = globalDefault;
+ }
+}
+
+void MVAvarsBase::SetWP(double WP){
+ btagMcut = WP;
+}
+
+void MVAvarsBase::FillMVAvarMap(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ std::cout << "WARNING: This is the base of 'FillMVAvarMap', this should never happen!\n";
+}
+
+void MVAvarsBase::FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ std::cout << "WARNING: This is the base of 'FillMVAvarMap', this should never happen!\n";
+}
+
+
+std::map<std::string, float> MVAvarsBase::GetMVAvars(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ FillMVAvarMap(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4);
+ return variableMap;
+}
+
+std::map<std::string, float> MVAvarsBase::GetMVAvars( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ FillMVAvarMap(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4, jets_idx);
+ return variableMap;
+}
+
+std::map<std::string, TLorentzVector> MVAvarsBase::GetVectors( const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ return std::map<std::string, TLorentzVector>();
+}
+
+bool MVAvarsBase::SkipEvent(const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx)
+{
+ std::cout << "CAREFUL - THIS IS THE BASE 'SkipEvent' FUNCTION!\n";
+ return false;
+}
+
+
+std::map<std::string, float>& MVAvarsBase::GetVariables()
+{
+ return variableMap;
+}
+
+float MVAvarsBase::GetHT(const std::vector<TLorentzVector> &selectedJetP4)
+{
+ float ht=0;
+ for(unsigned int i=0; i<selectedJetP4.size(); i++)
+ {
+ ht += selectedJetP4[i].Pt();
+ }
+
+ return ht;
+}
+
+
+TLorentzVector MVAvarsBase::GetLeptonicW(const TLorentzVector &leptonP4, const TLorentzVector &metP4)
+{
+ TLorentzVector newMETP4;
+
+ // reconstruct leptonic W Chwalek method
+ float nu_e = std::sqrt(metP4.Px()*metP4.Px()+metP4.Py()*metP4.Py());
+ float nu_px = metP4.Px();
+ float nu_py = metP4.Py();
+
+ float mw = 80.43;
+
+ float mtsq= ((leptonP4.Pt()+nu_e)*(leptonP4.Pt()+nu_e) -
+ (leptonP4.Px()+nu_px)*(leptonP4.Px()+nu_px) -
+ (leptonP4.Py()+nu_py)*(leptonP4.Py()+nu_py));
+
+ float mt = (mtsq>0.0) ? std::sqrt(mtsq) : 0.0;
+
+ float A, B, Csq, C, S1, S2;
+ float scf(1.0);
+
+ if (mt < mw) {
+ A = mw*mw/2.0;
+ }
+ else {
+ A = mt*mt/2.0;
+ float k = nu_e*leptonP4.Pt() - nu_px*leptonP4.Px() - nu_py*leptonP4.Py();
+ if (k <0.0001) k = 0.0001;
+ scf = 0.5*(mw*mw)/k;
+ nu_px *= scf;
+ nu_py *= scf;
+ nu_e = sqrt(nu_px*nu_px + nu_py*nu_py);
+ }
+
+ B = nu_px*leptonP4.Px() + nu_py*leptonP4.Py();
+ Csq= 1.0 + nu_e*nu_e*(leptonP4.Pz()*leptonP4.Pz()-leptonP4.E()*leptonP4.E())/(A+B)/(A+B);
+ C = (Csq>0.0) ? std::sqrt(Csq) : 0.0;
+ S1 = (-(A+B)*leptonP4.Pz() + (A+B)*leptonP4.E()*C)/(leptonP4.Pz()*leptonP4.Pz()-leptonP4.E()*leptonP4.E());
+ S2 = (-(A+B)*leptonP4.Pz() - (A+B)*leptonP4.E()*C)/(leptonP4.Pz()*leptonP4.Pz()-leptonP4.E()*leptonP4.E());
+
+ float nu_pz = (std::abs(S1) < std::abs(S2)) ? S1 : S2;
+ nu_e = sqrt(metP4.Px()*metP4.Px()+metP4.Py()*metP4.Py()+nu_pz*nu_pz);
+ newMETP4.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
+
+
+ return leptonP4 + newMETP4;
+}
\ No newline at end of file
diff --git a/src/MVAvarsJABDTthq.cpp b/src/MVAvarsJABDTthq.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9ba29ec6c607f8afb9816fad75215cd5f69baf4f
--- /dev/null
+++ b/src/MVAvarsJABDTthq.cpp
@@ -0,0 +1,204 @@
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTthq.h"
+
+using namespace std;
+
+MVAvarsJABDTthq::MVAvarsJABDTthq()
+{
+ Recolabels = { "Reco_tHq_top_m"
+ ,"Reco_tHq_top_pt"
+ ,"Reco_tHq_top_phi"
+ ,"Reco_tHq_top_eta"
+ ,"Reco_tHq_top_h_dr"
+
+ ,"Reco_tHq_btop_m"
+ ,"Reco_tHq_btop_pt"
+ ,"Reco_tHq_btop_phi"
+ ,"Reco_tHq_btop_eta"
+ ,"Reco_tHq_btop_idx"
+ ,"Reco_tHq_btop_lep_dr"
+ ,"Reco_tHq_btop_w_dr"
+
+ ,"Reco_tHq_ljet_m"
+ ,"Reco_tHq_ljet_pt"
+ ,"Reco_tHq_ljet_phi"
+ ,"Reco_tHq_ljet_eta"
+ ,"Reco_tHq_ljet_idx"
+
+ ,"Reco_tHq_h_m"
+ ,"Reco_tHq_h_pt"
+ ,"Reco_tHq_h_phi"
+ ,"Reco_tHq_h_eta"
+ ,"Reco_tHq_h_dr"
+
+ ,"Reco_tHq_hdau_m1"
+ ,"Reco_tHq_hdau_pt1"
+ ,"Reco_tHq_hdau_phi1"
+ ,"Reco_tHq_hdau_eta1"
+ ,"Reco_tHq_hdau_idx1"
+
+ ,"Reco_tHq_hdau_m2"
+ ,"Reco_tHq_hdau_pt2"
+ ,"Reco_tHq_hdau_phi2"
+ ,"Reco_tHq_hdau_eta2"
+ ,"Reco_tHq_hdau_idx2"
+
+ ,"Reco_tHq_costhetastar"
+
+ ,"Reco_JABDT_tHq_log_top_m"
+ ,"Reco_JABDT_tHq_log_top_pt"
+ ,"Reco_JABDT_tHq_log_h_m"
+ ,"Reco_JABDT_tHq_log_h_pt"
+ ,"Reco_JABDT_tHq_log_ljet_pt"
+ ,"Reco_JABDT_tHq_abs_top_eta"
+ ,"Reco_JABDT_tHq_abs_btop_eta"
+ ,"Reco_JABDT_tHq_abs_h_eta"
+ ,"Reco_JABDT_tHq_abs_ljet_eta"
+ ,"Reco_JABDT_tHq_abs_ljet_eta__M__btop_eta"
+ ,"Reco_JABDT_tHq_abs_top_eta__M__higg_eta"
+ ,"Reco_JABDT_tHq_Jet_CSV_hdau1"
+ ,"Reco_JABDT_tHq_Jet_CSV_hdau2"
+ ,"Reco_JABDT_tHq_Jet_CSV_btop"
+ ,"Reco_JABDT_tHq_Jet_CSV_ljet"
+ ,"Reco_JABDT_tHq_costheta_btop_lep"
+ ,"Reco_JABDT_tHq_log_min_hdau1_pt_hdau2_pt"
+ ,"Reco_JABDT_tHq_ljet_e__M__btop_e"
+ ,"Reco_JABDT_tHq_top_pt__P__h_pt__P__ljet_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"
+ };
+
+ for(uint i=0; i<Recolabels.size(); i++)
+ {
+ variableMap[Recolabels[i]] = globalDefault;
+ }
+}
+
+MVAvarsJABDTthq::~MVAvarsJABDTthq()
+{
+}
+
+void MVAvarsJABDTthq::FillMVAvarMap(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ // Reset all map entries to globalDefault value so that noting is left over from the last event
+ ResetVariableMap();
+
+ // ==================================================
+ // construct object vectors etc
+ float HT = GetHT(selectedJetP4);
+
+ std::map<std::string, TLorentzVector> vectors = GetVectors(selectedLeptonP4[0], selectedJetP4, metP4, jets_idx);
+
+
+ variableMap["Reco_tHq_btop_idx"] = jets_idx.at(tHqIndexes::tHq_btop_idx);
+ variableMap["Reco_tHq_ljet_idx"] = jets_idx.at(tHqIndexes::tHq_ljet_idx);
+ variableMap["Reco_tHq_hdau_idx1"] = jets_idx.at(tHqIndexes::tHq_hdau1_idx);
+ variableMap["Reco_tHq_hdau_idx2"] = jets_idx.at(tHqIndexes::tHq_hdau2_idx);
+
+ variableMap["Reco_tHq_top_m"] = vectors["top"].M();
+ variableMap["Reco_tHq_top_pt"] = vectors["top"].Pt();
+ variableMap["Reco_tHq_top_phi"] = vectors["top"].Phi();
+ variableMap["Reco_tHq_top_eta"] = vectors["top"].Eta();
+ variableMap["Reco_tHq_top_h_dr"] = vectors["top"].DeltaR(vectors["higg"]);
+
+ variableMap["Reco_tHq_btop_m"] = vectors["btop"].M();
+ variableMap["Reco_tHq_btop_pt"] = vectors["btop"].Pt();
+ variableMap["Reco_tHq_btop_phi"] = vectors["btop"].Phi();
+ variableMap["Reco_tHq_btop_eta"] = vectors["btop"].Eta();
+ variableMap["Reco_tHq_btop_lep_dr"]= vectors["btop"].DeltaR(selectedLeptonP4[0]);
+ variableMap["Reco_tHq_btop_w_dr"] = vectors["btop"].DeltaR(vectors["leptonicW"]);
+
+ variableMap["Reco_tHq_ljet_m"] = vectors["ljet"].M();
+ variableMap["Reco_tHq_ljet_pt"] = vectors["ljet"].Pt();
+ variableMap["Reco_tHq_ljet_phi"] = vectors["ljet"].Phi();
+ variableMap["Reco_tHq_ljet_eta"] = vectors["ljet"].Eta();
+
+ variableMap["Reco_tHq_h_m"] = vectors["higg"].M();
+ variableMap["Reco_tHq_h_pt"] = vectors["higg"].Pt();
+ variableMap["Reco_tHq_h_phi"] = vectors["higg"].Phi();
+ variableMap["Reco_tHq_h_eta"] = vectors["higg"].Eta();
+ variableMap["Reco_tHq_h_dr"] = vectors["hdau1"].DeltaR(vectors["hdau2"]);
+
+ variableMap["Reco_tHq_hdau_m1"] = vectors["hdau1"].M();
+ variableMap["Reco_tHq_hdau_pt1"] = vectors["hdau1"].Pt();
+ variableMap["Reco_tHq_hdau_phi1"] = vectors["hdau1"].Phi();
+ variableMap["Reco_tHq_hdau_eta1"] = vectors["hdau1"].Eta();
+
+ variableMap["Reco_tHq_hdau_m2"] = vectors["hdau2"].M();
+ variableMap["Reco_tHq_hdau_pt2"] = vectors["hdau2"].Pt();
+ variableMap["Reco_tHq_hdau_phi2"] = vectors["hdau2"].Phi();
+ variableMap["Reco_tHq_hdau_eta2"] = vectors["hdau2"].Eta();
+
+ //variables requiring a boost
+ TVector3 top_boost = vectors["top"].BoostVector();
+
+ TLorentzVector lep_p4_boost = selectedLeptonP4[0];
+ lep_p4_boost.Boost(-top_boost);
+
+ TLorentzVector ljet_boost = vectors["ljet"];
+ ljet_boost.Boost(-top_boost);
+ variableMap["Reco_tHq_costhetastar"] = lep_p4_boost.Vect().Dot(ljet_boost.Vect()) / (lep_p4_boost.Vect().Mag()*ljet_boost.Vect().Mag());
+
+ // variables for JABDT
+ variableMap["Reco_JABDT_tHq_log_top_m"] = log(vectors["top"].M());
+ variableMap["Reco_JABDT_tHq_log_top_pt"] = log(vectors["top"].Pt());
+ variableMap["Reco_JABDT_tHq_log_h_m"] = log(vectors["higg"].M());
+ variableMap["Reco_JABDT_tHq_log_h_pt"] = log(vectors["higg"].Pt());
+ variableMap["Reco_JABDT_tHq_log_ljet_pt"] = log(vectors["ljet"].Pt());
+
+ variableMap["Reco_JABDT_tHq_abs_top_eta"] = fabs(vectors["top"].Eta());
+ variableMap["Reco_JABDT_tHq_abs_btop_eta"] = fabs(vectors["btop"].Eta());
+ variableMap["Reco_JABDT_tHq_abs_h_eta"] = fabs(vectors["higg"].Eta());
+ variableMap["Reco_JABDT_tHq_abs_ljet_eta"] = fabs(vectors["ljet"].Eta());
+
+ variableMap["Reco_JABDT_tHq_abs_ljet_eta__M__btop_eta"] = fabs(vectors["ljet"].Eta()-vectors["btop"].Eta());
+ variableMap["Reco_JABDT_tHq_abs_top_eta__M__higg_eta"] = fabs(vectors["top"].Eta()-vectors["higg"].Eta());
+
+ variableMap["Reco_JABDT_tHq_Jet_CSV_hdau1"] = selectedJetCSV[jets_idx.at(tHqIndexes::tHq_hdau1_idx)];
+ variableMap["Reco_JABDT_tHq_Jet_CSV_hdau2"] = selectedJetCSV[jets_idx.at(tHqIndexes::tHq_hdau2_idx)];
+ variableMap["Reco_JABDT_tHq_Jet_CSV_btop"] = selectedJetCSV[jets_idx.at(tHqIndexes::tHq_btop_idx)];
+ variableMap["Reco_JABDT_tHq_Jet_CSV_ljet"] = selectedJetCSV[jets_idx.at(tHqIndexes::tHq_ljet_idx)];
+
+ variableMap["Reco_JABDT_tHq_costheta_btop_lep"] = selectedLeptonP4[0].Vect().Dot(vectors["btop"].Vect()) / (selectedLeptonP4[0].Vect().Mag()*vectors["btop"].Vect().Mag());
+ variableMap["Reco_JABDT_tHq_log_min_hdau1_pt_hdau2_pt"] = log(min(vectors["hdau1"].Pt(), vectors["hdau2"].Pt()));
+ variableMap["Reco_JABDT_tHq_ljet_e__M__btop_e"] = vectors["ljet"].E() - vectors["btop"].E();
+ variableMap["Reco_JABDT_tHq_top_pt__P__h_pt__P__ljet_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"] = (vectors["top"].Pt() + vectors["higg"].Pt() + vectors["ljet"].Pt())/(HT + metP4.Pt() + selectedLeptonP4[0].Pt());
+}
+
+std::map<std::string, TLorentzVector> MVAvarsJABDTthq::GetVectors(const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ std::map<std::string, TLorentzVector> vectors;
+
+ vectors["leptonicW"] = GetLeptonicW(selectedLeptonP4, metP4);
+
+ vectors["btop"] = selectedJetP4[jets_idx.at(tHqIndexes::tHq_btop_idx)];
+ vectors["ljet"] = selectedJetP4[jets_idx.at(tHqIndexes::tHq_ljet_idx)];
+ vectors["hdau1"] = selectedJetP4[jets_idx.at(tHqIndexes::tHq_hdau1_idx)];
+ vectors["hdau2"] = selectedJetP4[jets_idx.at(tHqIndexes::tHq_hdau2_idx)];
+
+ vectors["top"] = vectors["btop"] + vectors["leptonicW"];
+ vectors["higg"] = vectors["hdau1"] + vectors["hdau2"];
+
+ return vectors;
+}
+
+bool MVAvarsJABDTthq::SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx)
+{
+ int btop = jets_idx.at(tHqIndexes::tHq_btop_idx);
+ int ljet = jets_idx.at(tHqIndexes::tHq_ljet_idx);
+ int hdau1 = jets_idx.at(tHqIndexes::tHq_hdau1_idx);
+ int hdau2 = jets_idx.at(tHqIndexes::tHq_hdau2_idx);
+
+ if (!JetIsCentral(selectedJetP4[btop])) return true;
+ if (!JetIsCentral(selectedJetP4[hdau1])) return true;
+ if (!JetIsCentral(selectedJetP4[hdau2])) return true;
+ if (JetIsTagged(selectedJetP4[ljet], selectedJetCSV[ljet]) || !JetIsSelected(selectedJetP4[ljet])) return true;
+
+ return false;
+}
\ No newline at end of file
diff --git a/src/MVAvarsJABDTthw.cpp b/src/MVAvarsJABDTthw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..555da5ffe513b18e63b96dbb4d3d31a0fa7424cf
--- /dev/null
+++ b/src/MVAvarsJABDTthw.cpp
@@ -0,0 +1,282 @@
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTthw.h"
+
+using namespace std;
+
+MVAvarsJABDTthw::MVAvarsJABDTthw()
+{
+ topisleptonic = false;
+
+ Recolabels = { "Reco_tHW_btop_idx"
+ ,"Reco_tHW_whaddau_idx1"
+ ,"Reco_tHW_whaddau_idx2"
+ ,"Reco_tHW_hdau_idx1"
+ ,"Reco_tHW_hdau_idx2"
+ ,"Reco_tHW_leptonictop"
+
+ ,"Reco_tHW_top_m"
+ ,"Reco_tHW_top_pt"
+ ,"Reco_tHW_top_phi"
+ ,"Reco_tHW_top_eta"
+ ,"Reco_tHW_top_h_dr"
+
+ ,"Reco_tHW_btop_m"
+ ,"Reco_tHW_btop_pt"
+ ,"Reco_tHW_btop_phi"
+ ,"Reco_tHW_btop_eta"
+ ,"Reco_tHW_btop_lepw_dr"
+
+ ,"Reco_tHW_whad_m"
+ ,"Reco_tHW_whad_pt"
+ ,"Reco_tHW_whad_phi"
+ ,"Reco_tHW_whad_eta"
+ ,"Reco_tHW_whad_dr"
+
+ ,"Reco_tHW_whaddau_m1"
+ ,"Reco_tHW_whaddau_pt1"
+ ,"Reco_tHW_whaddau_phi1"
+ ,"Reco_tHW_whaddau_eta1"
+
+ ,"Reco_tHW_whaddau_m2"
+ ,"Reco_tHW_whaddau_pt2"
+ ,"Reco_tHW_whaddau_phi2"
+ ,"Reco_tHW_whaddau_eta2"
+
+ ,"Reco_tHW_h_m"
+ ,"Reco_tHW_h_pt"
+ ,"Reco_tHW_h_phi"
+ ,"Reco_tHW_h_eta"
+ ,"Reco_tHW_h_dr"
+
+ ,"Reco_tHW_hdau_m1"
+ ,"Reco_tHW_hdau_pt1"
+ ,"Reco_tHW_hdau_phi1"
+ ,"Reco_tHW_hdau_eta1"
+
+ ,"Reco_tHW_hdau_m2"
+ ,"Reco_tHW_hdau_pt2"
+ ,"Reco_tHW_hdau_phi2"
+ ,"Reco_tHW_hdau_eta2"
+
+ ,"Reco_tHW_top_m"
+ ,"Reco_tHW_top_pt"
+ ,"Reco_tHW_top_phi"
+ ,"Reco_tHW_top_eta"
+ ,"Reco_tHW_top_h_dr"
+
+ ,"Reco_tHW_wtop_m"
+ ,"Reco_tHW_wtop_pt"
+ ,"Reco_tHW_wtop_phi"
+ ,"Reco_tHW_wtop_eta"
+ ,"Reco_tHW_wtop_h_dr"
+
+ ,"Reco_tHW_wb_m"
+ ,"Reco_tHW_wb_pt"
+ ,"Reco_tHW_wb_phi"
+ ,"Reco_tHW_wb_eta"
+ ,"Reco_tHW_wb_h_dr"
+
+ ,"Reco_JABDT_tHW_Jet_CSV_btop"
+ ,"Reco_JABDT_tHW_Jet_CSV_hdau1"
+ ,"Reco_JABDT_tHW_Jet_CSV_hdau2"
+
+ ,"Reco_JABDT_tHW_log_top_m"
+ ,"Reco_JABDT_tHW_log_top_pt"
+ ,"Reco_JABDT_tHW_log_h_m"
+ ,"Reco_JABDT_tHW_log_h_pt"
+ ,"Reco_JABDT_tHW_log_wb_m"
+ ,"Reco_JABDT_tHW_log_wb_pt"
+ ,"Reco_JABDT_tHW_log_whad_m"
+ ,"Reco_JABDT_tHW_log_whad_pt"
+ ,"Reco_JABDT_tHW_wlep_pt__M__whad_pt"
+ ,"Reco_JABDT_tHW_abs_wlep_eta__M__whad_eta"
+
+ ,"Reco_JABDT_tHW_abs_top_eta__M__wb_eta"
+ ,"Reco_JABDT_tHW_abs_btop_eta"
+ ,"Reco_JABDT_tHW_abs_top_eta__M__higg_eta"
+ ,"Reco_JABDT_tHW_abs_wb_eta"
+ ,"Reco_JABDT_tHW_abs_top_eta"
+ ,"Reco_JABDT_tHW_top_pt__P__h_pt__P__wb_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"
+ ,"Reco_JABDT_tHW_costheta_btop_lep"
+ };
+
+ for(uint i=0; i<Recolabels.size(); i++)
+ {
+ variableMap[Recolabels[i]] = globalDefault;
+ }
+}
+
+MVAvarsJABDTthw::~MVAvarsJABDTthw()
+{
+}
+
+void MVAvarsJABDTthw::FillMVAvarMap(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ // Reset all map entries to globalDefault value so that noting is left over from the last event
+ ResetVariableMap();
+
+ // ==================================================
+ // construct object vectors etc
+ float HT = GetHT(selectedJetP4);
+
+ std::map<std::string, TLorentzVector> vectors = GetVectors(selectedLeptonP4[0], selectedJetP4, metP4, jets_idx);
+
+ variableMap["Reco_tHW_btop_idx"] = jets_idx.at(tHWIndexes::tHW_btop_idx);
+ variableMap["Reco_tHW_whaddau_idx1"]= jets_idx.at(tHWIndexes::tHW_whaddau1_idx);
+ variableMap["Reco_tHW_whaddau_idx2"]= jets_idx.at(tHWIndexes::tHW_whaddau2_idx);
+ variableMap["Reco_tHW_hdau_idx1"] = jets_idx.at(tHWIndexes::tHW_hdau1_idx);
+ variableMap["Reco_tHW_hdau_idx2"] = jets_idx.at(tHWIndexes::tHW_hdau2_idx);
+ variableMap["Reco_tHW_leptonictop"] = topisleptonic;
+
+ variableMap["Reco_tHW_top_m"] = vectors["top"].M();
+ variableMap["Reco_tHW_top_pt"] = vectors["top"].Pt();
+ variableMap["Reco_tHW_top_phi"] = vectors["top"].Phi();
+ variableMap["Reco_tHW_top_eta"] = vectors["top"].Eta();
+ variableMap["Reco_tHW_top_h_dr"] = vectors["top"].DeltaR(vectors["higg"]);
+
+ variableMap["Reco_tHW_btop_m"] = vectors["btop"].M();
+ variableMap["Reco_tHW_btop_pt"] = vectors["btop"].Pt();
+ variableMap["Reco_tHW_btop_phi"] = vectors["btop"].Phi();
+ variableMap["Reco_tHW_btop_eta"] = vectors["btop"].Eta();
+ variableMap["Reco_tHW_btop_lepw_dr"]= vectors["btop"].DeltaR(vectors["leptonicW"]);
+
+ variableMap["Reco_tHW_whad_m"] = vectors["whad"].M();
+ variableMap["Reco_tHW_whad_pt"] = vectors["whad"].Pt();
+ variableMap["Reco_tHW_whad_phi"] = vectors["whad"].Phi();
+ variableMap["Reco_tHW_whad_eta"] = vectors["whad"].Eta();
+ variableMap["Reco_tHW_whad_dr"] = vectors["whaddau1"].DeltaR(vectors["whaddau2"]);
+
+ variableMap["Reco_tHW_whaddau_m1"] = vectors["whaddau1"].M();
+ variableMap["Reco_tHW_whaddau_pt1"] = vectors["whaddau1"].Pt();
+ variableMap["Reco_tHW_whaddau_phi1"]= vectors["whaddau1"].Phi();
+ variableMap["Reco_tHW_whaddau_eta1"]= vectors["whaddau1"].Eta();
+
+ variableMap["Reco_tHW_whaddau_m2"] = vectors["whaddau2"].M();
+ variableMap["Reco_tHW_whaddau_pt2"] = vectors["whaddau2"].Pt();
+ variableMap["Reco_tHW_whaddau_phi2"]= vectors["whaddau2"].Phi();
+ variableMap["Reco_tHW_whaddau_eta2"]= vectors["whaddau2"].Eta();
+
+ variableMap["Reco_tHW_h_m"] = vectors["higg"].M();
+ variableMap["Reco_tHW_h_pt"] = vectors["higg"].Pt();
+ variableMap["Reco_tHW_h_phi"] = vectors["higg"].Phi();
+ variableMap["Reco_tHW_h_eta"] = vectors["higg"].Eta();
+ variableMap["Reco_tHW_h_dr"] = vectors["hdau1"].DeltaR(vectors["hdau2"]);
+
+ variableMap["Reco_tHW_hdau_m1"] = vectors["hdau1"].M();
+ variableMap["Reco_tHW_hdau_pt1"] = vectors["hdau1"].Pt();
+ variableMap["Reco_tHW_hdau_phi1"] = vectors["hdau1"].Phi();
+ variableMap["Reco_tHW_hdau_eta1"] = vectors["hdau1"].Eta();
+
+ variableMap["Reco_tHW_hdau_m2"] = vectors["hdau2"].M();
+ variableMap["Reco_tHW_hdau_pt2"] = vectors["hdau2"].Pt();
+ variableMap["Reco_tHW_hdau_phi2"] = vectors["hdau2"].Phi();
+ variableMap["Reco_tHW_hdau_eta2"] = vectors["hdau2"].Eta();
+
+ variableMap["Reco_tHW_top_m"] = vectors["top"].M();
+ variableMap["Reco_tHW_top_pt"] = vectors["top"].Pt();
+ variableMap["Reco_tHW_top_phi"] = vectors["top"].Phi();
+ variableMap["Reco_tHW_top_eta"] = vectors["top"].Eta();
+ variableMap["Reco_tHW_top_h_dr"] = vectors["top"].DeltaR(vectors["higg"]);
+
+ variableMap["Reco_tHW_wtop_m"] = vectors["wtop"].M();
+ variableMap["Reco_tHW_wtop_pt"] = vectors["wtop"].Pt();
+ variableMap["Reco_tHW_wtop_phi"] = vectors["wtop"].Phi();
+ variableMap["Reco_tHW_wtop_eta"] = vectors["wtop"].Eta();
+ variableMap["Reco_tHW_wtop_h_dr"] = vectors["wtop"].DeltaR(vectors["higg"]);
+
+ variableMap["Reco_tHW_wb_m"] = vectors["wb"].M();
+ variableMap["Reco_tHW_wb_pt"] = vectors["wb"].Pt();
+ variableMap["Reco_tHW_wb_phi"] = vectors["wb"].Phi();
+ variableMap["Reco_tHW_wb_eta"] = vectors["wb"].Eta();
+ variableMap["Reco_tHW_wb_h_dr"] = vectors["wb"].DeltaR(vectors["higg"]);
+
+ // variables for JABDT
+ variableMap["Reco_JABDT_tHW_Jet_CSV_btop"] = selectedJetCSV[jets_idx.at(tHWIndexes::tHW_btop_idx)];
+ variableMap["Reco_JABDT_tHW_Jet_CSV_hdau1"] = selectedJetCSV[jets_idx.at(tHWIndexes::tHW_hdau1_idx)];
+ variableMap["Reco_JABDT_tHW_Jet_CSV_hdau2"] = selectedJetCSV[jets_idx.at(tHWIndexes::tHW_hdau2_idx)];
+
+ variableMap["Reco_JABDT_tHW_log_top_m"] = log(vectors["top"].M());
+ variableMap["Reco_JABDT_tHW_log_top_pt"] = log(vectors["top"].Pt());
+ variableMap["Reco_JABDT_tHW_log_h_m"] = log(vectors["higg"].M());
+ variableMap["Reco_JABDT_tHW_log_h_pt"] = log(vectors["higg"].Pt());
+ variableMap["Reco_JABDT_tHW_log_wb_m"] = log(vectors["wb"].M());
+ variableMap["Reco_JABDT_tHW_log_wb_pt"] = log(vectors["wb"].Pt());
+ variableMap["Reco_JABDT_tHW_log_whad_m"] = log(vectors["whad"].M());
+ variableMap["Reco_JABDT_tHW_log_whad_pt"] = log(vectors["whad"].Pt());
+ variableMap["Reco_JABDT_tHW_wlep_pt__M__whad_pt"] = vectors["leptonicW"].Pt() - vectors["whad"].Pt();
+ variableMap["Reco_JABDT_tHW_abs_wlep_eta__M__whad_eta"] = fabs(vectors["leptonicW"].Eta() - vectors["whad"].Eta());
+
+ variableMap["Reco_JABDT_tHW_abs_top_eta__M__wb_eta"] = fabs(vectors["top"].Eta() - vectors["wb"].Eta());
+ variableMap["Reco_JABDT_tHW_abs_btop_eta"] = fabs(vectors["btop"].Eta());
+ variableMap["Reco_JABDT_tHW_abs_top_eta__M__higg_eta"] = fabs(vectors["top"].Eta() - vectors["higg"].Eta());
+ variableMap["Reco_JABDT_tHW_abs_wb_eta"] = fabs(vectors["wb"].Eta());
+ variableMap["Reco_JABDT_tHW_abs_top_eta"] = fabs(vectors["top"].Eta());
+
+ variableMap["Reco_JABDT_tHW_top_pt__P__h_pt__P__wb_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"] = (vectors["top"].Pt() + vectors["higg"].Pt() + vectors["wb"].Pt())/(HT + metP4.Pt() + selectedLeptonP4[0].Pt());
+ variableMap["Reco_JABDT_tHW_costheta_btop_lep"] = selectedLeptonP4[0].Vect().Dot(vectors["btop"].Vect()) / (selectedLeptonP4[0].Vect().Mag()*vectors["btop"].Vect().Mag());
+}
+
+std::map<std::string, TLorentzVector> MVAvarsJABDTthw::GetVectors(const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ std::map<std::string, TLorentzVector> vectors;
+
+ vectors["leptonicW"] = GetLeptonicW(selectedLeptonP4, metP4);
+
+ vectors["btop"] = selectedJetP4[jets_idx.at(tHWIndexes::tHW_btop_idx)];
+ vectors["whaddau1"] = selectedJetP4[jets_idx.at(tHWIndexes::tHW_whaddau1_idx)];
+ vectors["whaddau2"] = selectedJetP4[jets_idx.at(tHWIndexes::tHW_whaddau2_idx)];
+ vectors["hdau1"] = selectedJetP4[jets_idx.at(tHWIndexes::tHW_hdau1_idx)];
+ vectors["hdau2"] = selectedJetP4[jets_idx.at(tHWIndexes::tHW_hdau2_idx)];
+
+ vectors["higg"] = vectors["hdau1"] + vectors["hdau2"];
+ vectors["whad"] = vectors["whaddau1"] + vectors["whaddau2"];
+
+ vectors["toplep"] = vectors["leptonicW"] + vectors["btop"];
+ vectors["tophad"] = vectors["whad"] + vectors["btop"];
+
+ if(topisleptonic)
+ {
+ vectors["wtop"] = vectors["leptonicW"];
+ vectors["wb"] = vectors["whad"];
+ }
+ else
+ {
+ vectors["wtop"] = vectors["whad"];
+ vectors["wb"] = vectors["leptonicW"];
+ }
+ vectors["top"] = vectors["wtop"] + vectors["btop"];
+
+ return vectors;
+}
+
+bool MVAvarsJABDTthw::SkipEvent(const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx)
+{
+ int btop = jets_idx.at(tHWIndexes::tHW_btop_idx);
+ int whaddau1 = jets_idx.at(tHWIndexes::tHW_whaddau1_idx);
+ int whaddau2 = jets_idx.at(tHWIndexes::tHW_whaddau2_idx);
+ int hdau1 = jets_idx.at(tHWIndexes::tHW_hdau1_idx);
+ int hdau2 = jets_idx.at(tHWIndexes::tHW_hdau2_idx);
+
+ if (!JetIsCentral(selectedJetP4[btop])) return true;
+ if (!JetIsCentral(selectedJetP4[hdau1])) return true;
+ if (!JetIsCentral(selectedJetP4[hdau2])) return true;
+ if (JetIsTagged(selectedJetP4[whaddau1], selectedJetCSV.at(whaddau1))) return true;
+ if (JetIsTagged(selectedJetP4[whaddau2], selectedJetCSV.at(whaddau2))) return true;
+
+ int btags = 0;
+ if(JetIsTagged(selectedJetP4.at(btop), selectedJetCSV.at(btop))) btags++;
+ if(JetIsTagged(selectedJetP4.at(hdau1), selectedJetCSV.at(hdau1))) btags++;
+ if(JetIsTagged(selectedJetP4.at(hdau2), selectedJetCSV.at(hdau2))) btags++;
+
+ if(btags < 1) return true;
+
+ return false;
+}
\ No newline at end of file
diff --git a/src/MVAvarsJABDTttbar.cpp b/src/MVAvarsJABDTttbar.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cc2b1bc8ac670c155eb17d3bcdf9c0cd7e843f89
--- /dev/null
+++ b/src/MVAvarsJABDTttbar.cpp
@@ -0,0 +1,192 @@
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTttbar.h"
+
+using namespace std;
+
+MVAvarsJABDTttbar::MVAvarsJABDTttbar()
+{
+ // ==================================================
+ //init all variables used for ttbar hypothesis testing
+ Recolabels = { "Reco_ttbar_btoplep_m"
+ ,"Reco_ttbar_btoplep_pt"
+ ,"Reco_ttbar_btoplep_phi"
+ ,"Reco_ttbar_btoplep_eta"
+ ,"Reco_ttbar_btoplep_idx"
+
+ ,"Reco_ttbar_whad_m"
+ ,"Reco_ttbar_whad_pt"
+ ,"Reco_ttbar_whad_phi"
+ ,"Reco_ttbar_whad_eta"
+ ,"Reco_ttbar_whad_dr"
+
+ ,"Reco_ttbar_whaddau_m1"
+ ,"Reco_ttbar_whaddau_m2"
+ ,"Reco_ttbar_whaddau_pt1"
+ ,"Reco_ttbar_whaddau_pt2"
+ ,"Reco_ttbar_whaddau_phi1"
+ ,"Reco_ttbar_whaddau_phi2"
+
+ ,"Reco_ttbar_whaddau_eta1"
+ ,"Reco_ttbar_whaddau_eta2"
+ ,"Reco_ttbar_whaddau_idx1"
+ ,"Reco_ttbar_whaddau_idx2"
+
+ ,"Reco_ttbar_btophad_m"
+ ,"Reco_ttbar_btophad_pt"
+ ,"Reco_ttbar_btophad_phi"
+ ,"Reco_ttbar_btophad_eta"
+ ,"Reco_ttbar_btophad_idx"
+
+ ,"Reco_ttbar_tophad_m"
+ ,"Reco_ttbar_tophad_pt"
+ ,"Reco_ttbar_tophad_phi"
+ ,"Reco_ttbar_tophad_eta"
+ ,"Reco_ttbar_tophad_dr"
+
+ ,"Reco_ttbar_toplep_m"
+ ,"Reco_ttbar_toplep_pt"
+ ,"Reco_ttbar_toplep_phi"
+ ,"Reco_ttbar_toplep_eta"
+
+ ,"Reco_LeptonicW_Eta"
+ ,"Reco_LeptonicW_M"
+ ,"Reco_LeptonicW_Phi"
+ ,"Reco_LeptonicW_Pt"
+
+ ,"Reco_JABDT_ttbar_log_whad_m"
+ ,"Reco_JABDT_ttbar_Jet_CSV_btoplep"
+ ,"Reco_JABDT_ttbar_Jet_CSV_btophad"
+ ,"Reco_JABDT_ttbar_Jet_CSV_whaddau1"
+ ,"Reco_JABDT_ttbar_Jet_CSV_whaddau2"
+ ,"Reco_JABDT_ttbar_log_tophad_m__M__whad_m"
+ ,"Reco_JABDT_ttbar_log_tophad_pt"
+ ,"Reco_JABDT_ttbar_log_toplep_m"
+ ,"Reco_JABDT_ttbar_log_toplep_pt"
+ ,"Reco_JABDT_ttbar_tophad_pt__P__toplep_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"
+ ,"Reco_JABDT_ttbar_costheta_toplep_tophad"
+ };
+
+ for(uint i=0; i<Recolabels.size(); i++)
+ {
+ variableMap[Recolabels[i]] = globalDefault;
+ }
+}
+
+MVAvarsJABDTttbar::~MVAvarsJABDTttbar()
+{
+}
+
+void MVAvarsJABDTttbar::FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ // Reset all map entries to globalDefault value so that noting is left over from the last event
+ ResetVariableMap();
+ // ==================================================
+ // construct object vectors etc
+ float HT = GetHT(selectedJetP4);
+
+ std::map<std::string, TLorentzVector> vectors = GetVectors(selectedLeptonP4[0], selectedJetP4, metP4, jets_idx);
+
+ variableMap["Reco_ttbar_btoplep_idx"] = jets_idx.at(ttbarIndexes::ttbar_btoplep_idx);
+ variableMap["Reco_ttbar_btophad_idx"] = jets_idx.at(ttbarIndexes::ttbar_btophad_idx);
+ variableMap["Reco_ttbar_whaddau_idx1"] = jets_idx.at(ttbarIndexes::ttbar_whaddau1_idx);
+ variableMap["Reco_ttbar_whaddau_idx2"] = jets_idx.at(ttbarIndexes::ttbar_whaddau2_idx);
+
+ variableMap["Reco_ttbar_btoplep_m"] = vectors["btoplep"].M();
+ variableMap["Reco_ttbar_btoplep_pt"] = vectors["btoplep"].Pt();
+ variableMap["Reco_ttbar_btoplep_phi"] = vectors["btoplep"].Phi();
+ variableMap["Reco_ttbar_btoplep_eta"] = vectors["btoplep"].Eta();
+
+ variableMap["Reco_ttbar_whad_m"] = vectors["whad"].M();
+ variableMap["Reco_ttbar_whad_pt"] = vectors["whad"].Pt();
+ variableMap["Reco_ttbar_whad_phi"] = vectors["whad"].Phi();
+ variableMap["Reco_ttbar_whad_eta"] = vectors["whad"].Eta();
+ variableMap["Reco_ttbar_whad_dr"] = vectors["whaddau1"].DeltaR(vectors["whaddau2"]);
+
+ variableMap["Reco_ttbar_whaddau_m1"] = vectors["whaddau1"].M();
+ variableMap["Reco_ttbar_whaddau_pt1"] = vectors["whaddau1"].Pt();
+ variableMap["Reco_ttbar_whaddau_phi1"] = vectors["whaddau1"].Phi();
+ variableMap["Reco_ttbar_whaddau_eta1"] = vectors["whaddau1"].Eta();
+
+ variableMap["Reco_ttbar_whaddau_m2"] = vectors["whaddau2"].M();
+ variableMap["Reco_ttbar_whaddau_pt2"] = vectors["whaddau2"].Pt();
+ variableMap["Reco_ttbar_whaddau_phi2"] = vectors["whaddau2"].Phi();
+ variableMap["Reco_ttbar_whaddau_eta2"] = vectors["whaddau2"].Eta();
+
+ variableMap["Reco_ttbar_btophad_m"] = vectors["btophad"].M();
+ variableMap["Reco_ttbar_btophad_pt"] = vectors["btophad"].Pt();
+ variableMap["Reco_ttbar_btophad_phi"] = vectors["btophad"].Phi();
+ variableMap["Reco_ttbar_btophad_eta"] = vectors["btophad"].Eta();
+
+ variableMap["Reco_ttbar_tophad_m"] = vectors["tophad"].M();
+ variableMap["Reco_ttbar_tophad_pt"] = vectors["tophad"].Pt();
+ variableMap["Reco_ttbar_tophad_phi"] = vectors["tophad"].Phi();
+ variableMap["Reco_ttbar_tophad_eta"] = vectors["tophad"].Eta();
+ variableMap["Reco_ttbar_tophad_dr"] = vectors["whad"].DeltaR(vectors["btophad"]);
+
+ variableMap["Reco_ttbar_toplep_m"] = vectors["toplep"].M();
+ variableMap["Reco_ttbar_toplep_pt"] = vectors["toplep"].Pt();
+ variableMap["Reco_ttbar_toplep_phi"] = vectors["toplep"].Phi();
+ variableMap["Reco_ttbar_toplep_eta"] = vectors["toplep"].Eta();
+
+ variableMap["Reco_LeptonicW_Eta"] = vectors["leptonicW"].Eta();
+ variableMap["Reco_LeptonicW_M"] = vectors["leptonicW"].M();
+ variableMap["Reco_LeptonicW_Phi"] = vectors["leptonicW"].Phi();
+ variableMap["Reco_LeptonicW_Pt"] = vectors["leptonicW"].Pt();
+
+ variableMap["Reco_JABDT_ttbar_log_whad_m"] = log(vectors["whad"].M());
+
+ variableMap["Reco_JABDT_ttbar_Jet_CSV_btoplep"] = selectedJetCSV[jets_idx.at(ttbarIndexes::ttbar_btoplep_idx)];
+ variableMap["Reco_JABDT_ttbar_Jet_CSV_btophad"] = selectedJetCSV[jets_idx.at(ttbarIndexes::ttbar_btophad_idx)];
+ variableMap["Reco_JABDT_ttbar_Jet_CSV_whaddau1"] = selectedJetCSV[jets_idx.at(ttbarIndexes::ttbar_whaddau1_idx)];
+ variableMap["Reco_JABDT_ttbar_Jet_CSV_whaddau2"] = selectedJetCSV[jets_idx.at(ttbarIndexes::ttbar_whaddau2_idx)];
+
+ variableMap["Reco_JABDT_ttbar_log_tophad_m__M__whad_m"] = log(vectors["tophad"].M() - vectors["whad"].M());
+ variableMap["Reco_JABDT_ttbar_log_tophad_pt"] = log(vectors["tophad"].Pt());
+
+ variableMap["Reco_JABDT_ttbar_log_toplep_m"] = log(vectors["toplep"].M());
+ variableMap["Reco_JABDT_ttbar_log_toplep_pt"] = log(vectors["toplep"].Pt());
+
+ variableMap["Reco_JABDT_ttbar_tophad_pt__P__toplep_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"] = (vectors["tophad"].Pt() + vectors["toplep"].Pt())/(HT + metP4.Pt() + selectedLeptonP4[0].Pt());
+ variableMap["Reco_JABDT_ttbar_costheta_toplep_tophad"] = vectors["toplep"].Vect().Dot(vectors["tophad"].Vect()) / (vectors["toplep"].Vect().Mag()*vectors["tophad"].Vect().Mag());
+}
+
+std::map<std::string, TLorentzVector> MVAvarsJABDTttbar::GetVectors(const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ std::map<std::string, TLorentzVector> vectors;
+
+ vectors["leptonicW"] = GetLeptonicW(selectedLeptonP4, metP4);
+
+ vectors["btoplep"] = selectedJetP4[jets_idx.at(ttbarIndexes::ttbar_btoplep_idx)];
+ vectors["btophad"] = selectedJetP4[jets_idx.at(ttbarIndexes::ttbar_btophad_idx)];
+ vectors["whaddau1"] = selectedJetP4[jets_idx.at(ttbarIndexes::ttbar_whaddau1_idx)];
+ vectors["whaddau2"] = selectedJetP4[jets_idx.at(ttbarIndexes::ttbar_whaddau2_idx)];
+
+ vectors["toplep"] = vectors["btoplep"] + vectors["leptonicW"];
+ vectors["whad"] = vectors["whaddau1"] + vectors["whaddau2"];
+ vectors["tophad"] = vectors["whad"] + vectors["btophad"];
+
+ return vectors;
+}
+
+bool MVAvarsJABDTttbar::SkipEvent( const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx)
+{
+ int btoplep = jets_idx.at(ttbarIndexes::ttbar_btoplep_idx);
+ int btophad = jets_idx.at(ttbarIndexes::ttbar_btophad_idx);
+ int whaddau1 = jets_idx.at(ttbarIndexes::ttbar_whaddau1_idx);
+ int whaddau2 = jets_idx.at(ttbarIndexes::ttbar_whaddau2_idx);
+
+ if (!JetIsTagged(selectedJetP4[btophad], selectedJetCSV[btophad])) return true;
+ if (!JetIsTagged(selectedJetP4[btoplep], selectedJetCSV[btoplep])) return true;
+ if (!JetIsSelected(selectedJetP4[whaddau1])) return true;
+ if (!JetIsSelected(selectedJetP4[whaddau2])) return true;
+
+ return false;
+}
\ No newline at end of file
diff --git a/src/MVAvarsJABDTtth.cpp b/src/MVAvarsJABDTtth.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..330d4782b035f896d54525bb023360c2c7d8c0a1
--- /dev/null
+++ b/src/MVAvarsJABDTtth.cpp
@@ -0,0 +1,240 @@
+#include "TTH/CommonClassifier/interface/MVAvarsJABDTtth.h"
+
+using namespace std;
+
+MVAvarsJABDTtth::MVAvarsJABDTtth()
+{
+ // ==================================================
+ //init all variables used for ttH hypothesis testing
+ Recolabels = { "Reco_ttH_btoplep_m"
+ ,"Reco_ttH_btoplep_pt"
+ ,"Reco_ttH_btoplep_phi"
+ ,"Reco_ttH_btoplep_eta"
+ ,"Reco_ttH_btoplep_idx"
+ ,"Reco_ttH_btoplep_w_dr"
+
+ ,"Reco_ttH_whad_m"
+ ,"Reco_ttH_whad_pt"
+ ,"Reco_ttH_whad_phi"
+ ,"Reco_ttH_whad_eta"
+ ,"Reco_ttH_whad_dr"
+
+ ,"Reco_ttH_whaddau_m1"
+ ,"Reco_ttH_whaddau_pt1"
+ ,"Reco_ttH_whaddau_phi1"
+ ,"Reco_ttH_whaddau_eta1"
+ ,"Reco_ttH_whaddau_idx1"
+
+ ,"Reco_ttH_whaddau_m2"
+ ,"Reco_ttH_whaddau_pt2"
+ ,"Reco_ttH_whaddau_phi2"
+ ,"Reco_ttH_whaddau_eta2"
+ ,"Reco_ttH_whaddau_idx2"
+
+ ,"Reco_ttH_btophad_m"
+ ,"Reco_ttH_btophad_pt"
+ ,"Reco_ttH_btophad_phi"
+ ,"Reco_ttH_btophad_eta"
+ ,"Reco_ttH_btophad_idx"
+
+ ,"Reco_ttH_tophad_m"
+ ,"Reco_ttH_tophad_pt"
+ ,"Reco_ttH_tophad_phi"
+ ,"Reco_ttH_tophad_eta"
+ ,"Reco_ttH_tophad_dr"
+
+ ,"Reco_ttH_toplep_m"
+ ,"Reco_ttH_toplep_pt"
+ ,"Reco_ttH_toplep_phi"
+ ,"Reco_ttH_toplep_eta"
+ ,"Reco_ttH_toplep_w_dr"
+
+ ,"Reco_ttH_h_m"
+ ,"Reco_ttH_h_pt"
+ ,"Reco_ttH_h_phi"
+ ,"Reco_ttH_h_eta"
+ ,"Reco_ttH_h_dr"
+
+ ,"Reco_ttH_hdau_m1"
+ ,"Reco_ttH_hdau_pt1"
+ ,"Reco_ttH_hdau_phi1"
+ ,"Reco_ttH_hdau_eta1"
+ ,"Reco_ttH_hdau_m2"
+
+ ,"Reco_ttH_hdau_pt2"
+ ,"Reco_ttH_hdau_phi2"
+ ,"Reco_ttH_hdau_eta2"
+ ,"Reco_ttH_hdau_idx1"
+ ,"Reco_ttH_hdau_idx2"
+
+ ,"Reco_JABDT_ttH_log_whad_m"
+ ,"Reco_JABDT_ttH_Jet_CSV_btophad"
+ ,"Reco_JABDT_ttH_Jet_CSV_btoplep"
+ ,"Reco_JABDT_ttH_Jet_CSV_hdau1"
+ ,"Reco_JABDT_ttH_Jet_CSV_hdau2"
+ ,"Reco_JABDT_ttH_Jet_CSV_whaddau1"
+ ,"Reco_JABDT_ttH_Jet_CSV_whaddau2"
+ ,"Reco_JABDT_ttH_log_tophad_m__M__whad_m"
+ ,"Reco_JABDT_ttH_log_tophad_pt"
+ ,"Reco_JABDT_ttH_log_toplep_m"
+ ,"Reco_JABDT_ttH_log_toplep_pt"
+ ,"Reco_JABDT_ttH_tophad_pt__P__toplep_pt__P__h_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"
+ ,"Reco_JABDT_ttH_log_h_pt"
+ ,"Reco_JABDT_ttH_log_h_m"
+ };
+
+ for(uint i=0; i<Recolabels.size(); i++)
+ {
+ variableMap[Recolabels[i]] = globalDefault;
+ }
+}
+
+MVAvarsJABDTtth::~MVAvarsJABDTtth()
+{
+}
+
+void MVAvarsJABDTtth::FillMVAvarMap( const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ // Reset all map entries to globalDefault value so that noting is left over from the last event
+ ResetVariableMap();
+
+ // ==================================================
+ // construct object vectors etc
+ float HT = GetHT(selectedJetP4);
+ std::map<std::string, TLorentzVector> vectors = GetVectors(selectedLeptonP4[0], selectedJetP4, metP4, jets_idx);
+
+ variableMap["Reco_ttH_btoplep_idx"] = jets_idx.at(ttHIndexes::ttH_btoplep_idx);
+ variableMap["Reco_ttH_btophad_idx"] = jets_idx.at(ttHIndexes::ttH_btophad_idx);
+ variableMap["Reco_ttH_whaddau_idx1"] = jets_idx.at(ttHIndexes::ttH_whaddau1_idx);
+ variableMap["Reco_ttH_whaddau_idx2"] = jets_idx.at(ttHIndexes::ttH_whaddau2_idx);
+ variableMap["Reco_ttH_hdau_idx1"] = jets_idx.at(ttHIndexes::ttH_hdau1_idx);
+ variableMap["Reco_ttH_hdau_idx2"] = jets_idx.at(ttHIndexes::ttH_hdau2_idx);
+
+ variableMap["Reco_ttH_btoplep_m"] = vectors["btoplep"].M();
+ variableMap["Reco_ttH_btoplep_pt"] = vectors["btoplep"].Pt();
+ variableMap["Reco_ttH_btoplep_phi"] = vectors["btoplep"].Phi();
+ variableMap["Reco_ttH_btoplep_eta"] = vectors["btoplep"].Eta();
+ variableMap["Reco_ttH_btoplep_w_dr"] = vectors["btoplep"].DeltaR(vectors["leptonicWP4"]);
+
+ variableMap["Reco_ttH_whad_m"] = vectors["whad"].M();
+ variableMap["Reco_ttH_whad_pt"] = vectors["whad"].Pt();
+ variableMap["Reco_ttH_whad_phi"] = vectors["whad"].Phi();
+ variableMap["Reco_ttH_whad_eta"] = vectors["whad"].Eta();
+ variableMap["Reco_ttH_whad_dr"] = vectors["whaddau1"].DeltaR(vectors["whaddau2"]);
+
+ variableMap["Reco_ttH_whaddau_m1"] = vectors["whaddau1"].M();
+ variableMap["Reco_ttH_whaddau_pt1"] = vectors["whaddau1"].Pt();
+ variableMap["Reco_ttH_whaddau_phi1"] = vectors["whaddau1"].Phi();
+ variableMap["Reco_ttH_whaddau_eta1"] = vectors["whaddau1"].Eta();
+
+ variableMap["Reco_ttH_whaddau_m2"] = vectors["whaddau2"].M();
+ variableMap["Reco_ttH_whaddau_pt2"] = vectors["whaddau2"].Pt();
+ variableMap["Reco_ttH_whaddau_phi2"] = vectors["whaddau2"].Phi();
+ variableMap["Reco_ttH_whaddau_eta2"] = vectors["whaddau2"].Eta();
+
+ variableMap["Reco_ttH_btophad_m"] = vectors["btophad"].M();
+ variableMap["Reco_ttH_btophad_pt"] = vectors["btophad"].Pt();
+ variableMap["Reco_ttH_btophad_phi"] = vectors["btophad"].Phi();
+ variableMap["Reco_ttH_btophad_eta"] = vectors["btophad"].Eta();
+
+ variableMap["Reco_ttH_tophad_m"] = vectors["tophad"].M();
+ variableMap["Reco_ttH_tophad_pt"] = vectors["tophad"].Pt();
+ variableMap["Reco_ttH_tophad_phi"] = vectors["tophad"].Phi();
+ variableMap["Reco_ttH_tophad_eta"] = vectors["tophad"].Eta();
+ variableMap["Reco_ttH_tophad_dr"] = vectors["whad"].DeltaR(vectors["btophad"]);
+
+ variableMap["Reco_ttH_toplep_m"] = vectors["toplep"].M();
+ variableMap["Reco_ttH_toplep_pt"] = vectors["toplep"].Pt();
+ variableMap["Reco_ttH_toplep_phi"] = vectors["toplep"].Phi();
+ variableMap["Reco_ttH_toplep_eta"] = vectors["toplep"].Eta();
+ variableMap["Reco_ttH_toplep_w_dr"] = vectors["toplep"].DeltaR(vectors["leptonicWP4"]);
+
+ variableMap["Reco_ttH_h_m"] = vectors["higg"].M();
+ variableMap["Reco_ttH_h_pt"] = vectors["higg"].Pt();
+ variableMap["Reco_ttH_h_phi"] = vectors["higg"].Phi();
+ variableMap["Reco_ttH_h_eta"] = vectors["higg"].Eta();
+ variableMap["Reco_ttH_h_dr"] = vectors["hdau1"].DeltaR(vectors["hdau2"]);
+
+ variableMap["Reco_ttH_hdau_m1"] = vectors["hdau1"].M();
+ variableMap["Reco_ttH_hdau_pt1"] = vectors["hdau1"].Pt();
+ variableMap["Reco_ttH_hdau_phi1"] = vectors["hdau1"].Phi();
+ variableMap["Reco_ttH_hdau_eta1"] = vectors["hdau1"].Eta();
+
+ variableMap["Reco_ttH_hdau_m2"] = vectors["hdau2"].M();
+ variableMap["Reco_ttH_hdau_pt2"] = vectors["hdau2"].Pt();
+ variableMap["Reco_ttH_hdau_phi2"] = vectors["hdau2"].Phi();
+ variableMap["Reco_ttH_hdau_eta2"] = vectors["hdau2"].Eta();
+
+ // variables for JABDT
+ variableMap["Reco_JABDT_ttH_log_whad_m"] = log(vectors["whad"].M());
+ variableMap["Reco_JABDT_ttH_Jet_CSV_btoplep"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_btoplep_idx)];
+ variableMap["Reco_JABDT_ttH_Jet_CSV_btophad"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_btophad_idx)];
+ variableMap["Reco_JABDT_ttH_Jet_CSV_whaddau1"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_whaddau1_idx)];
+ variableMap["Reco_JABDT_ttH_Jet_CSV_whaddau2"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_whaddau2_idx)];
+ variableMap["Reco_JABDT_ttH_Jet_CSV_hdau1"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_hdau1_idx)];
+ variableMap["Reco_JABDT_ttH_Jet_CSV_hdau2"] = selectedJetCSV[jets_idx.at(ttHIndexes::ttH_hdau2_idx)];
+ variableMap["Reco_JABDT_ttH_log_tophad_m__M__whad_m"] = log(vectors["tophad"].M() - vectors["whad"].M());
+ variableMap["Reco_JABDT_ttH_log_tophad_pt"] = log(vectors["tophad"].Pt());
+ variableMap["Reco_JABDT_ttH_log_toplep_m"] = log(vectors["toplep"].M());
+ variableMap["Reco_JABDT_ttH_log_toplep_pt"] = log(vectors["toplep"].Pt());
+ variableMap["Reco_JABDT_ttH_tophad_pt__P__toplep_pt__P__h_pt__DIV__Evt_HT__P__Evt_Pt_MET__P__Lep_Pt"] = (vectors["tophad"].Pt() + vectors["toplep"].Pt() + vectors["higg"].Pt())/(HT + metP4.Pt() + selectedLeptonP4[0].Pt());
+ variableMap["Reco_JABDT_ttH_log_h_pt"] = log(vectors["higg"].Pt());
+ variableMap["Reco_JABDT_ttH_log_h_m"] = log(vectors["higg"].M());
+}
+
+std::map<std::string, TLorentzVector> MVAvarsJABDTtth::GetVectors(const TLorentzVector &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const TLorentzVector &metP4,
+ const std::vector<int> &jets_idx)
+{
+ std::map<std::string, TLorentzVector> vectors;
+
+ vectors["leptonicW"] = GetLeptonicW(selectedLeptonP4, metP4);
+
+ vectors["btoplep"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_btoplep_idx)];
+ vectors["btophad"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_btophad_idx)];
+ vectors["whaddau1"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_whaddau1_idx)];
+ vectors["whaddau2"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_whaddau2_idx)];
+ vectors["hdau1"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_hdau1_idx)];
+ vectors["hdau2"] = selectedJetP4[jets_idx.at(ttHIndexes::ttH_hdau2_idx)];
+
+ vectors["toplep"] = vectors["btoplep"] + vectors["leptonicW"];
+ vectors["whad"] = vectors["whaddau1"] + vectors["whaddau2"];
+ vectors["tophad"] = vectors["whad"] + vectors["btophad"];
+ vectors["higg"] = vectors["hdau1"] + vectors["hdau2"];
+
+ return vectors;
+}
+
+bool MVAvarsJABDTtth::SkipEvent(const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const std::vector<int> &jets_idx)
+{
+ int btoplep = jets_idx.at(ttHIndexes::ttH_btoplep_idx);
+ int btophad = jets_idx.at(ttHIndexes::ttH_btophad_idx);
+ int hdau1 = jets_idx.at(ttHIndexes::ttH_hdau1_idx);
+ int hdau2 = jets_idx.at(ttHIndexes::ttH_hdau2_idx);
+ int whaddau1 = jets_idx.at(ttHIndexes::ttH_whaddau1_idx);
+ int whaddau2 = jets_idx.at(ttHIndexes::ttH_whaddau2_idx);
+
+ if (!JetIsCentral(selectedJetP4.at(btoplep))) return true;
+ if (!JetIsCentral(selectedJetP4.at(btophad))) return true;
+ if (!JetIsCentral(selectedJetP4.at(hdau1))) return true;
+ if (!JetIsCentral(selectedJetP4.at(hdau2))) return true;
+ if (!JetIsSelected(selectedJetP4.at(whaddau1))) return true;
+ if (!JetIsSelected(selectedJetP4.at(whaddau2))) return true;
+
+ int btags = 0;
+ if(JetIsTagged(selectedJetP4.at(btophad), selectedJetCSV.at(btophad))) btags++;
+ if(JetIsTagged(selectedJetP4.at(btoplep), selectedJetCSV.at(btoplep))) btags++;
+ if(JetIsTagged(selectedJetP4.at(hdau1), selectedJetCSV.at(hdau1))) btags++;
+ if(JetIsTagged(selectedJetP4.at(hdau2), selectedJetCSV.at(hdau2))) btags++;
+
+ if(btags < 2) return true;
+
+ return false;
+}
\ No newline at end of file
diff --git a/src/PermutationManager.cpp b/src/PermutationManager.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c2548462a8e8f97d70c18a4d43b792f007f2823d
--- /dev/null
+++ b/src/PermutationManager.cpp
@@ -0,0 +1,156 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// PermutationManager
+// ---------------
+//
+// 13/05/2019 Marco Link
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include "../interface/PermutationManager.h"
+
+
+// #include <cassert>
+#include <iostream>
+#include <fstream>
+// #include <sstream>
+
+
+using namespace std;
+
+
+////////////////////////////////////////////////////////////////////////////////
+// construction / destruction
+////////////////////////////////////////////////////////////////////////////////
+
+//______________________________________________________________________________
+PermutationManager::PermutationManager(unsigned int N_idx, unsigned int Jets_Max) : Nidx(N_idx), JetsMax(Jets_Max)
+{
+ if(Nidx > JetsMax)
+ {
+ cout << "error: number of indizes > max number of jets" << endl;
+ }
+
+ for(unsigned int Njets=Nidx; Njets<=JetsMax; Njets++)
+ {
+ std::vector<std::vector<int>> temp;
+
+ unsigned int comb[Nidx];
+ for(unsigned int idx=0; idx<Nidx; idx++)
+ {
+ comb[idx] = Njets-1;
+ }
+
+ unsigned int sum = 999;
+ bool reduce;
+ bool skip;
+
+ // fill permutations
+ while(sum > 0)
+ {
+ std::vector<int> combination;
+ reduce = true;
+ skip = false;
+ sum = 0;
+ for(unsigned int idx=0; idx<Nidx; idx++)
+ {
+ combination.push_back(comb[idx]);
+
+ // skip combinations with double indizes
+ for(unsigned int idx2=0; idx2<idx; idx2++)
+ {
+ if(comb[idx] == comb[idx2]) skip=true;
+ }
+ }
+ if(!skip) temp.push_back(combination);
+
+ for(unsigned int idx=0; idx<Nidx; idx++)
+ {
+ // reduce by 1
+ if(reduce && comb[idx]>0)
+ {
+ comb[idx] = comb[idx] - 1;
+ reduce = false;
+ }
+ else if (reduce)
+ {
+ comb[idx] = Njets-1;
+ }
+
+ sum += comb[idx];
+ }
+ }
+ permutations.push_back(temp);
+ }
+}
+
+
+//______________________________________________________________________________
+PermutationManager::~PermutationManager()
+{
+
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// implementation of member functions
+////////////////////////////////////////////////////////////////////////////////
+unsigned int PermutationManager::get_Npermutations(unsigned int Njets)
+{
+ return permutations.at(getIndex(Njets)).size();
+}
+
+
+void PermutationManager::get_permutation(std::vector<int>* permutation, unsigned int Njets, unsigned int i)
+{
+ if(i >= get_Npermutations(Njets)) cout << "error: permutation index out of range!" << endl;
+
+ *permutation = permutations[getIndex(Njets)].at(i);
+}
+
+
+void PermutationManager::show()
+{
+ cout << "Permutationmanager for " << Nidx << " indizes" << endl;
+ cout << "NJets (Index), theo, generated permutations" << endl;
+ for(unsigned int i=Nidx; i<=JetsMax; i++)
+ {
+ cout << setw(2) << i << "(" << setw(2) << getIndex(i) << "), " << setw(15) << TMath::Factorial(i)/TMath::Factorial(i - Nidx) << ", " << setw(15) << get_Npermutations(i) << endl;
+ }
+ cout << endl;
+}
+
+
+void PermutationManager::show(unsigned int Njets)
+{
+ std::vector<std::vector<int>> combinations = permutations.at(getIndex(Njets));
+
+ cout << "permutations for " << Nidx << " indizes and " << Njets << " jets" << endl;
+
+ for(unsigned int n=0; n<get_Npermutations(Njets); n++)
+ {
+ for(unsigned int idx=0; idx<Nidx; idx++)
+ {
+ cout << combinations.at(n).at(idx) << "\t";
+ }
+ cout << endl;
+ }
+ cout << endl;
+}
+
+
+unsigned int PermutationManager::getIndex(unsigned int Njets)
+{
+ if(Njets > JetsMax)
+ {
+ cout << "error: permutation not generated for " << Njets << " increase maximum jet number!" << endl;
+ return 0;
+ }
+ else if(Njets < Nidx)
+ {
+ cout << "error: jet number " << Njets << " smaller than number of indizes! Automagically fixing problem ... NOT" << endl;
+ return 0;
+ }
+
+ return Njets - Nidx;
+}
diff --git a/src/ReconstructedVars.cpp b/src/ReconstructedVars.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f68a44ff8c6bda820f52b42d52edb41f207074e
--- /dev/null
+++ b/src/ReconstructedVars.cpp
@@ -0,0 +1,161 @@
+#include "TTH/CommonClassifier/interface/ReconstructedVars.h"
+
+using namespace std;
+
+ReconstructedVars::ReconstructedVars(bool ttH, bool ttZ):
+ reconstruct_ttH(ttH), reconstruct_ttZ(ttZ),
+ chi2reco_ttH("RecoTTH_", "Higgs", 116.3, 17.6),
+ chi2reco_ttZ("RecoTTZ_", "Z", 86.8, 13.5)
+{
+ // ==================================================
+ //init all variables potentially used
+ if(reconstruct_ttH)
+ {
+ // angles of top system
+ std::map<std::string, double> topSystemAngles = chi2reco_ttH.GetVariableMap_TopSystemAngles();
+ for (auto it = topSystemAngles.begin(); it != topSystemAngles.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // objects of top system
+ std::map<std::string, double> topSystemObjects = chi2reco_ttH.GetVariableMap_TopSystemObjects();
+ for (auto it = topSystemObjects.begin(); it != topSystemObjects.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // angles of higgs system
+ std::map<std::string, double> higgsSystemAngles = chi2reco_ttH.GetVariableMap_BosonSystemAngles();
+ for (auto it = higgsSystemAngles.begin(); it != higgsSystemAngles.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // objects of higgs system
+ std::map<std::string, double> higgsSystemObjects = chi2reco_ttH.GetVariableMap_BosonSystemObjects();
+ for (auto it = higgsSystemObjects.begin(); it != higgsSystemObjects.end(); it++)
+ variableMap[it->first] = -999.;
+ }
+
+
+ if(reconstruct_ttZ)
+ {
+ // angles of top system
+ std::map<std::string, double> topSystemAngles = chi2reco_ttZ.GetVariableMap_TopSystemAngles();
+ for (auto it = topSystemAngles.begin(); it != topSystemAngles.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // objects of top system
+ std::map<std::string, double> topSystemObjects = chi2reco_ttZ.GetVariableMap_TopSystemObjects();
+ for (auto it = topSystemObjects.begin(); it != topSystemObjects.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // angles of higgs system
+ std::map<std::string, double> ZSystemAngles = chi2reco_ttZ.GetVariableMap_BosonSystemAngles();
+ for (auto it = ZSystemAngles.begin(); it != ZSystemAngles.end(); it++)
+ variableMap[it->first] = -999.;
+
+ // objects of higgs system
+ std::map<std::string, double> ZSystemObjects = chi2reco_ttZ.GetVariableMap_BosonSystemObjects();
+ for (auto it = ZSystemObjects.begin(); it != ZSystemObjects.end(); it++)
+ variableMap[it->first] = -999.;
+ }
+}
+
+ReconstructedVars::~ReconstructedVars() {}
+
+
+void ReconstructedVars::ResetVariableMap()
+{
+ for (auto it = variableMap.begin(); it != variableMap.end(); it++)
+ it->second = -999.;
+}
+
+void ReconstructedVars::SetWP(double WP){
+ btagMcut = WP;
+}
+
+
+
+
+std::map<std::string, double> ReconstructedVars::GetReconstructedVars(
+ const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ //check if CSVWP is set properly
+ if(btagMcut<0){
+ std::cerr << "Please set CSV Wp properly\n" << std::endl;
+ throw std::exception();
+ }
+ // Reset all map entries to -999 so that noting is left over from the last event
+ ResetVariableMap();
+
+ // get number of jets
+ int nJets = selectedJetP4.size();
+
+ // Start Reconstruction
+ if(reconstruct_ttH)
+ {
+ chi2reco_ttH.ReconstructTTXSystem(selectedLeptonP4[0], selectedJetP4, selectedJetCSV, metP4, btagMcut);
+
+ ///////////// fill map /////////////
+ // angles of top system
+ std::map<std::string, double> topSystemAngles = chi2reco_ttH.FillVariableMap_TopSystemAngles();
+ for (auto it = topSystemAngles.begin(); it != topSystemAngles.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // objects of top system
+ std::map<std::string, double> topSystemObjects = chi2reco_ttH.FillVariableMap_TopSystemObjects();
+ for (auto it = topSystemObjects.begin(); it != topSystemObjects.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // only fill higgs system variables if more/equal 6 jets
+ if (nJets >= 6)
+ {
+ // angles of higgs system
+ std::map<std::string, double> higgsSystemAngles = chi2reco_ttH.FillVariableMap_BosonSystemAngles();
+ for (auto it = higgsSystemAngles.begin(); it != higgsSystemAngles.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // objects of higgs system
+ std::map<std::string, double> higgsSystemObjects = chi2reco_ttH.FillVariableMap_BosonSystemObjects();
+ for (auto it = higgsSystemObjects.begin(); it != higgsSystemObjects.end(); it++)
+ variableMap[it->first] = it->second;
+ }
+ }
+
+ if(reconstruct_ttZ)
+ {
+ chi2reco_ttZ.ReconstructTTXSystem(selectedLeptonP4[0], selectedJetP4, selectedJetCSV, metP4, btagMcut);
+
+ ///////////// fill map /////////////
+ // angles of top system
+ std::map<std::string, double> topSystemAngles = chi2reco_ttZ.FillVariableMap_TopSystemAngles();
+ for (auto it = topSystemAngles.begin(); it != topSystemAngles.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // objects of top system
+ std::map<std::string, double> topSystemObjects = chi2reco_ttZ.FillVariableMap_TopSystemObjects();
+ for (auto it = topSystemObjects.begin(); it != topSystemObjects.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // only fill higgs system variables if more/equal 6 jets
+ if (nJets >= 6)
+ {
+ // angles of higgs system
+ std::map<std::string, double> ZSystemAngles = chi2reco_ttZ.FillVariableMap_BosonSystemAngles();
+ for (auto it = ZSystemAngles.begin(); it != ZSystemAngles.end(); it++)
+ variableMap[it->first] = it->second;
+
+ // objects of higgs system
+ std::map<std::string, double> ZSystemObjects = chi2reco_ttZ.FillVariableMap_BosonSystemObjects();
+ for (auto it = ZSystemObjects.begin(); it != ZSystemObjects.end(); it++)
+ variableMap[it->first] = it->second;
+ }
+ }
+ return variableMap;
+}
+
+
+std::map<std::string, double> ReconstructedVars::GetVariables()
+{
+ ResetVariableMap();
+ return variableMap;
+}
diff --git a/src/thqHypothesisCombinatorics.cpp b/src/thqHypothesisCombinatorics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..047b6cd9c4fdef965d8b2fa5f4b19fcefd696ff5
--- /dev/null
+++ b/src/thqHypothesisCombinatorics.cpp
@@ -0,0 +1,26 @@
+#include "TTH/CommonClassifier/interface/thqHypothesisCombinatorics.h"
+
+
+thqHypothesisCombinatorics::thqHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring):
+HypothesisCombinatorics(optional_varstring)
+{
+ minJets = 4;
+ reader_th.reset(new TMVA::Reader( "!Color:!Silent" ));
+ permutator.reset(new PermutationManager(minJets, 10));
+ permutator->show();
+ mvars.reset(new MVAvarsJABDTthq());
+
+ bdt_name = "tHq JABDT";
+ bdtoutput_name = "Reco_tHq_bestJABDToutput";
+
+ // setup input variables for TMVA factory
+ FillVariableNameList(optional_varstring, BDTlabels);
+ for (unsigned ivar=0;ivar<BDTlabels.size();ivar++){
+ reader_th->AddVariable(BDTlabels.at(ivar).c_str(), mvars->GetAdress(BDTlabels[ivar]));
+ cout << "added variable to JABDT:" << BDTlabels.at(ivar).c_str() <<endl;
+ }
+
+ reader_th->BookMVA(bdt_name, weightpath);
+}
+
+thqHypothesisCombinatorics::~thqHypothesisCombinatorics(){}
\ No newline at end of file
diff --git a/src/thwHypothesisCombinatorics.cpp b/src/thwHypothesisCombinatorics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ea71b55aaf24faed7fd6e69b8ac0959e3d6fcb58
--- /dev/null
+++ b/src/thwHypothesisCombinatorics.cpp
@@ -0,0 +1,48 @@
+#include "TTH/CommonClassifier/interface/thwHypothesisCombinatorics.h"
+
+
+thwHypothesisCombinatorics::thwHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring):
+HypothesisCombinatorics(optional_varstring)
+{
+ minJets = 5;
+ reader_th.reset(new TMVA::Reader( "!Color:!Silent" ));
+ permutator.reset(new PermutationManager(minJets, 10));
+ permutator->show();
+ mvars.reset(new MVAvarsJABDTthw());
+
+ bdt_name = "tHW JABDT";
+ bdtoutput_name = "Reco_tHW_bestJABDToutput";
+
+ // setup input variables for TMVA factory
+ FillVariableNameList(optional_varstring, BDTlabels);
+ for (unsigned ivar=0;ivar<BDTlabels.size();ivar++){
+ reader_th->AddVariable(BDTlabels.at(ivar).c_str(), mvars->GetAdress(BDTlabels[ivar]));
+ cout << "added variable to JABDT:" << BDTlabels.at(ivar).c_str() <<endl;
+ }
+
+ reader_th->BookMVA(bdt_name, weightpath);
+}
+
+thwHypothesisCombinatorics::~thwHypothesisCombinatorics(){}
+
+std::map<std::string, float> thwHypothesisCombinatorics::GetBestPermutation(const std::vector<TLorentzVector> &selectedLeptonP4,
+ const std::vector<TLorentzVector> &selectedJetP4,
+ const std::vector<double> &selectedJetCSV,
+ const TLorentzVector &metP4)
+{
+ // reconstruct leptonic and hadronic top cases
+ mvars->SetTopisleptonic(true);
+ std::map<std::string, float> leptonictop = HypothesisCombinatorics::GetBestPermutation(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4);
+
+ mvars->SetTopisleptonic(false);
+ std::map<std::string, float> hadronictop = HypothesisCombinatorics::GetBestPermutation(selectedLeptonP4, selectedJetP4, selectedJetCSV, metP4);
+
+ if(hadronictop[bdtoutput_name] > leptonictop[bdtoutput_name])
+ {
+ return hadronictop;
+ }
+ else
+ {
+ return leptonictop;
+ }
+}
\ No newline at end of file
diff --git a/src/ttbarHypothesisCombinatorics.cpp b/src/ttbarHypothesisCombinatorics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4346a5a18ce084d775d30ae679b85540a316228c
--- /dev/null
+++ b/src/ttbarHypothesisCombinatorics.cpp
@@ -0,0 +1,26 @@
+#include "TTH/CommonClassifier/interface/ttbarHypothesisCombinatorics.h"
+
+
+ttbarHypothesisCombinatorics::ttbarHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring):
+HypothesisCombinatorics(optional_varstring)
+{
+ minJets = 4;
+ reader_th.reset(new TMVA::Reader( "!Color:!Silent" ));
+ permutator.reset(new PermutationManager(minJets, 10));
+ permutator->show();
+ mvars.reset(new MVAvarsJABDTttbar());
+
+ bdt_name = "ttbar JABDT";
+ bdtoutput_name = "Reco_ttbar_bestJABDToutput";
+
+ // setup input variables for TMVA factory
+ FillVariableNameList(optional_varstring, BDTlabels);
+ for (unsigned ivar=0;ivar<BDTlabels.size();ivar++){
+ reader_th->AddVariable(BDTlabels.at(ivar).c_str(), mvars->GetAdress(BDTlabels[ivar]));
+ cout << "added variable to JABDT:" << BDTlabels.at(ivar).c_str() <<endl;
+ }
+
+ reader_th->BookMVA(bdt_name, weightpath);
+}
+
+ttbarHypothesisCombinatorics::~ttbarHypothesisCombinatorics(){}
\ No newline at end of file
diff --git a/src/tthHypothesisCombinatorics.cpp b/src/tthHypothesisCombinatorics.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..141b00667f659c4ff1110492c595f49e04b14ed4
--- /dev/null
+++ b/src/tthHypothesisCombinatorics.cpp
@@ -0,0 +1,26 @@
+#include "TTH/CommonClassifier/interface/tthHypothesisCombinatorics.h"
+
+
+tthHypothesisCombinatorics::tthHypothesisCombinatorics(const std::string& weightpath, const std::string& optional_varstring):
+HypothesisCombinatorics(optional_varstring)
+{
+ minJets = 6;
+ reader_th.reset(new TMVA::Reader( "!Color:!Silent" ));
+ permutator.reset(new PermutationManager(minJets, 10));
+ permutator->show();
+ mvars.reset(new MVAvarsJABDTtth());
+
+ bdt_name = "ttH JABDT";
+ bdtoutput_name = "Reco_ttH_bestJABDToutput";
+
+ // setup input variables for TMVA factory
+ FillVariableNameList(optional_varstring, BDTlabels);
+ for (unsigned ivar=0;ivar<BDTlabels.size();ivar++){
+ reader_th->AddVariable(BDTlabels.at(ivar).c_str(), mvars->GetAdress(BDTlabels[ivar]));
+ cout << "added variable to JABDT:" << BDTlabels.at(ivar).c_str() <<endl;
+ }
+
+ reader_th->BookMVA(bdt_name, weightpath);
+}
+
+tthHypothesisCombinatorics::~tthHypothesisCombinatorics(){}
\ No newline at end of file