diff --git a/data/ci/template-analysis_shapeInterp_womcstats.txt b/data/ci/template-analysis_shapeInterp_womcstats.txt new file mode 100644 index 00000000000..cdb93252b8c --- /dev/null +++ b/data/ci/template-analysis_shapeInterp_womcstats.txt @@ -0,0 +1,89 @@ +# Datacard produced by CombineHarvester with git status: 8fe0e3c-dirty +imax 1 number of bins +jmax 6 number of processes minus 1 +kmax * number of nuisance parameters +-------------------------------------------------------------------------------- +shapes * htt_tt_9_13TeV htt_input.root htt_tt_9_13TeV/$PROCESS htt_tt_9_13TeV/$PROCESS_$SYSTEMATIC +shapes bbH htt_tt_9_13TeV htt_input.root htt_tt_9_13TeV/bbH$MASS htt_tt_9_13TeV/bbH$MASS_$SYSTEMATIC +shapes ggH htt_tt_9_13TeV htt_input.root htt_tt_9_13TeV/ggH$MASS htt_tt_9_13TeV/ggH$MASS_$SYSTEMATIC +-------------------------------------------------------------------------------- +bin htt_tt_9_13TeV +observation 3416.0 +-------------------------------------------------------------------------------- +bin htt_tt_9_13TeV htt_tt_9_13TeV htt_tt_9_13TeV htt_tt_9_13TeV htt_tt_9_13TeV htt_tt_9_13TeV htt_tt_9_13TeV +process ZL TTT VVT ZTT jetFakes ggH bbH +process 1 2 3 4 5 -1 0 +rate 37.5448 683.017 96.5185 742.649 2048.94 19.9504 198.521 +-------------------------------------------------------------------------------- +CMS_eff_b_13TeV lnN - 0.99/1.01 0.98/1.01 0.98/1.02 - 0.99/1.01 0.98/1.02 +CMS_eff_m lnN - - - 0.96 - - - +CMS_eff_t_13TeV lnN - 1.08 1.08 1.08 - 1.08 1.08 +CMS_eff_t_mssmHigh_tt_13TeV shape - 1 1 1 - 1 1 +CMS_eff_t_tt_13TeV lnN - 1.092 1.092 1.092 - 1.092 1.092 +CMS_fake_b_13TeV lnN 0.99/1.05 - 0.99/1.01 0.97/1.02 - 0.97/1.03 - +CMS_htt_dyShape_scale_m_13TeV shape - - - 1 - - - +CMS_htt_dyShape_stat_m400pt0_13TeV shape - - - 1 - - - +CMS_htt_dyShape_stat_m400pt40_13TeV shape - - - 1 - - - +CMS_htt_dyShape_stat_m400pt80_13TeV shape - - - 1 - - - +CMS_htt_dyShape_tjXsec_13TeV shape - - - 1 - - - +CMS_htt_eFakeTau_loose_13TeV lnN 1.03 - - - - - - +CMS_htt_mFakeTau_loose_13TeV lnN 1.05 - - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_10 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_11 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_12 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_13 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_14 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_15 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_16 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_8 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_TTT_bin_9 shape - 1 - - - - - +CMS_htt_tt_btag_13TeV_VVT_bin_17 shape - - 1 - - - - +CMS_htt_tt_btag_13TeV_VVT_bin_18 shape - - 1 - - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_1 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_15 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_17 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_18 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_2 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_3 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_4 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_ZTT_bin_5 shape - - - 1 - - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_1 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_10 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_11 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_12 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_13 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_14 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_16 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_17 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_2 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_3 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_5 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_6 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_7 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_8 shape - - - - 1 - - +CMS_htt_tt_btag_13TeV_jetFakes_bin_9 shape - - - - 1 - - +CMS_htt_ttbarAccept_tt_btag_13TeV lnN - 1.004 - - - - - +CMS_htt_vvXsec_13TeV lnN - - 1.05 - - - - +CMS_htt_zjXsec_13TeV lnN 1.04 - - - - - - +CMS_htt_zttAccept_tt_btag_13TeV lnN - - - 1.05 - - - +CMS_scale_j_13TeV lnN 0.96/1.01 - 0.99/1.01 0.99/1.01 - 0.99/1.01 0.99/1.01 +CMS_scale_t_1prong0pi0_13TeV shape - 1 1 1 - 1 1 +CMS_scale_t_1prong1pi0_13TeV shape - 1 1 1 - 1 1 +CMS_scale_t_3prong0pi0_13TeV shape - 1 1 1 - 1 1 +QCDScale_QshScale_bbH lnN - - - - - - 0.902 +ff_norm_stat_tt_tt_btag lnN - - - - 1.028 - - +ff_norm_syst_tt lnN - - - - 1.1 - - +ff_sub_syst_tt_tt_btag lnN - - - - 1.03 - - +lumi_13TeV lnN 1.025 - 1.025 - - 1.025 1.025 +norm_ff_dy_frac_tt_syst shape - - - - 1 - - +norm_ff_qcd_dm0_njet0_tt_stat shape - - - - 1 - - +norm_ff_qcd_dm0_njet1_tt_stat shape - - - - 1 - - +norm_ff_qcd_dm1_njet0_tt_stat shape - - - - 1 - - +norm_ff_qcd_dm1_njet1_tt_stat shape - - - - 1 - - +norm_ff_qcd_tt_syst shape - - - - 1 - - +norm_ff_tt_frac_tt_syst shape - - - - 1 - - +norm_ff_tt_tt_syst shape - - - - 1 - - +norm_ff_w_frac_tt_syst shape - - - - 1 - - +norm_ff_w_tt_syst shape - - - - 1 - - +rate_TT rateParam * TTT 1 [0,5] +rate_ZMM_ZTT_btag rateParam * ZTT 1.02 [0.8,1.2] diff --git a/docs/root_ad_development_notes.md b/docs/root_ad_development_notes.md new file mode 100644 index 00000000000..73c3a18ecfd --- /dev/null +++ b/docs/root_ad_development_notes.md @@ -0,0 +1,136 @@ +# Introduction +This file is created to keep track of issues and tests conducted for integrating AD in Combine. +## Setup + +Using a CMSSW build with ROOT 6.32.06: + +``` +cmssw-el8 +cmsrel CMSSW_14_2_0_pre2_ROOT632 +cd CMSSW_14_2_0_pre2_ROOT632/src +git clone https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit.git HiggsAnalysis/CombinedLimit +git checkout roofit_ad_ichep_2024_63206_comp +cd ../../ +cmsenv +scram b -j 8 +``` +## Tested models +Below you can find test results for several models with classes implemented in Combine. + +### Discrete profiling + +``` + text2workspace.py data/ci/datacard_RooMultiPdf.txt.gz -o multipdf_ws.root + python3 scripts/fitRooFitAD.py --input multipdf_ws.root +``` +Fails because only `RooAbsReal` parameters are allowed, pdf indices in `RooMultiPdf` models are `RooCategory`. +
+ Error message + +```bash + RooAbsReal* RooAbsPdf::createNLL(RooAbsData& data, const RooLinkedList& cmdArgs) => + runtime_error: In creation of function nll_func_wrapper wrapper: input param expected to be of type RooAbsReal. +``` + +
+ +### autoMCStats + +``` +text2workspace.py data/ci/template-analysis_shapeInterp.txt -o template_autoMCstats_ws.root -m 200 +python3 scripts/fitRooFitAD.py --input template_autoMCstats_ws.root +``` +Fails because AD is not implemented for "RooRealSumPdf" objects used to model MC statistical (https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/main/python/ShapeTools.py#L264) uncertainties with `autoMCstats` +
+ Error message + +```bash +[#0] ERROR:Minimization -- An analytical integral function for class "RooRealSumPdf" has not yet been implemented. +Traceback (most recent call last): + File "/afs/cern.ch/work/a/anigamov/CMSSW_14_2_ROOT632_X_2024-10-06-2300/src/HiggsAnalysis/CombinedLimit/scripts/fitRooFitAD.py", line 29, in + nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, EvalBackend="codegen") + File "/cvmfs/cms-ib.cern.ch/sw/x86_64/nweek-02858/el8_amd64_gcc12/lcg/root/6.32.06-f59fcaa47c786e7268e714e7e477ee41/lib/ROOT/_pythonization/_roofit/_rooabspdf.py", line 116, in createNLL + return self._createNLL["RooLinkedList const&"](args[0], _pack_cmd_args(*args[1:], **kwargs)) +cppyy.gbl.std.runtime_error: Could not find "createNLL" (set cppyy.set_debug() for C++ errors): + RooAbsReal* RooAbsPdf::createNLL(RooAbsData& data, const RooLinkedList& cmdArgs) => + runtime_error: An analytical integral function for class "RooRealSumPdf" has not yet been implemented. +``` +
+ +### Template based model without MC statistical uncertainties (autoMCStats) + +``` +text2workspace.py data/ci/template-analysis_shapeInterp_womcstats.txt -o template_ws.root -m 200 +python3 scripts/fitRooFitAD.py --input template_ws.root +``` +The fit runs with the warning shown below, but results are reasonable +
+ Warning: + +```bash +In module 'RooFitCore': +/cvmfs/cms-ib.cern.ch/sw/x86_64/nweek-02858/el8_amd64_gcc12/lcg/root/6.32.06-f59fcaa47c786e7268e714e7e477ee41/include/RooFit/Detail/MathFuncs.h:365:45: warning: function 'LnGamma' was not differentiated because clad failed to differentiate it and no suitable overload + was found in namespace 'custom_derivatives' + return pdf - weight * std::log(pdf) + TMath::LnGamma(weight + 1); + ^ +/cvmfs/cms-ib.cern.ch/sw/x86_64/nweek-02858/el8_amd64_gcc12/lcg/root/6.32.06-f59fcaa47c786e7268e714e7e477ee41/include/RooFit/Detail/MathFuncs.h:365:45: note: falling back to numerical differentiation for 'LnGamma' since no suitable overload was found and clad + could not derive it; to disable this feature, compile your programs with -DCLAD_NO_NUM_DIFF +``` + +
+ +### RooParametricHist + +``` +text2workspace.py data/ci/datacard_RooParametricHist.txt -o ws_RooParametricHist.root +python3 scripts/fitRooFitAD.py --input ws_RooParametricHist.root +``` + +
+ Error message + +```bash +[#0] ERROR:InputArguments -- RooHistPdf::weight(shapeBkg_tqq_muonCRfail2016) ERROR: Code Squashing currently only supports uniformly binned cases. +.... + ^ +[#0] ERROR:InputArguments -- Function roo_func_wrapper_0 could not be compiled. See above for details. +Traceback (most recent call last): + File "/afs/cern.ch/work/a/anigamov/CMSSW_14_2_ROOT632_X_2024-10-06-2300/src/HiggsAnalysis/CombinedLimit/scripts/fitRooFitAD.py", line 29, in + nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, EvalBackend="codegen") + File "/cvmfs/cms-ib.cern.ch/sw/x86_64/nweek-02858/el8_amd64_gcc12/lcg/root/6.32.06-f59fcaa47c786e7268e714e7e477ee41/lib/ROOT/_pythonization/_roofit/_rooabspdf.py", line 116, in createNLL + return self._createNLL["RooLinkedList const&"](args[0], _pack_cmd_args(*args[1:], **kwargs)) +cppyy.gbl.std.runtime_error: Could not find "createNLL" (set cppyy.set_debug() for C++ errors): + RooAbsReal* RooAbsPdf::createNLL(RooAbsData& data, const RooLinkedList& cmdArgs) => + runtime_error: Function roo_func_wrapper_0 could not be compiled. See above for details. +``` + +
+ +### RooHistPdf + +``` +ulimit -s unlimited +text2workspace.py data/ci/datacard_RooHistPdf.txt.gz -o ws_RooHistPdf.root +python3 scripts/fitRooFitAD.py --input ws_RooHistPdf.root +``` +AD is not implemented for `VerticalInterpPdf` class: https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/roofit_ad_ichep_2024_63206_comp/src/VerticalInterpPdf.cc +
+ Error message + +```bash +[#0] ERROR:Minimization -- An analytical integral function for class "VerticalInterpPdf" has not yet been implemented. +Traceback (most recent call last): + File "/afs/cern.ch/work/a/anigamov/CMSSW_14_2_ROOT632_X_2024-10-06-2300/src/HiggsAnalysis/CombinedLimit/scripts/fitRooFitAD.py", line 29, in + nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, EvalBackend="codegen") + File "/cvmfs/cms-ib.cern.ch/sw/x86_64/nweek-02858/el8_amd64_gcc12/lcg/root/6.32.06-f59fcaa47c786e7268e714e7e477ee41/lib/ROOT/_pythonization/_roofit/_rooabspdf.py", line 116, in createNLL + return self._createNLL["RooLinkedList const&"](args[0], _pack_cmd_args(*args[1:], **kwargs)) +cppyy.gbl.std.runtime_error: Could not find "createNLL" (set cppyy.set_debug() for C++ errors): + RooAbsReal* RooAbsPdf::createNLL(RooAbsData& data, const RooLinkedList& cmdArgs) => + runtime_error: An analytical integral function for class "VerticalInterpPdf" has not yet been implemented. +``` + +
+ +### Channel masking + +Channel masking is implemented [within eval of custom `CachingNLL` class](https://github.com/cms-analysis/HiggsAnalysis-CombinedLimit/blob/main/src/CachingNLL.cc#L1078). diff --git a/interface/AsymPow.h b/interface/AsymPow.h index 8346bdadba6..552c9829941 100644 --- a/interface/AsymPow.h +++ b/interface/AsymPow.h @@ -4,6 +4,7 @@ #include #include +#include "CombineCodegenImpl.h" //_________________________________________________ /* @@ -26,9 +27,13 @@ class AsymPow : public RooAbsReal { AsymPow(const char *name, const char *title, RooAbsReal &kappaLow, RooAbsReal &kappaHigh, RooAbsReal &theta) ; AsymPow(const AsymPow &other, const char *newname=0) ; - ~AsymPow() override ; + TObject * clone(const char *newname) const override { return new AsymPow(*this, newname); } - TObject * clone(const char *newname) const override ; + COMBINE_DECLARE_TRANSLATE; + + RooAbsReal const &kappaLow() const { return kappaLow_.arg(); } + RooAbsReal const &kappaHigh() const { return kappaHigh_.arg(); } + RooAbsReal const &theta() const { return theta_.arg(); } protected: Double_t evaluate() const override; @@ -40,4 +45,6 @@ class AsymPow : public RooAbsReal { ClassDefOverride(AsymPow,1) // Asymmetric power }; +COMBINE_DECLARE_CODEGEN_IMPL(AsymPow); + #endif diff --git a/interface/CombineCodegenImpl.h b/interface/CombineCodegenImpl.h new file mode 100644 index 00000000000..d1795eb1520 --- /dev/null +++ b/interface/CombineCodegenImpl.h @@ -0,0 +1,18 @@ +#ifndef HiggsAnalysis_CombinedLimit_CombineCodegenImpl_h +#define HiggsAnalysis_CombinedLimit_CombineCodegenImpl_h + +#include // for ROOT_VERSION + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +# define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME) \ + namespace RooFit { namespace Experimental { void codegenImpl(CLASS_NAME &arg, CodegenContext &ctx); }} +# define COMBINE_DECLARE_TRANSLATE +#elif ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +# define COMBINE_DECLARE_CODEGEN_IMPL(CLASS_NAME) +# define COMBINE_DECLARE_TRANSLATE void translate(RooFit::Detail::CodeSquashContext &ctx) const override; +#else +# define COMBINE_DECLARE_CODEGEN_IMPL(_) +# define COMBINE_DECLARE_TRANSLATE +#endif + +#endif diff --git a/interface/HZZ4LRooPdfs.h b/interface/HZZ4LRooPdfs.h index df1f7e2e3ae..fbae926a4d8 100644 --- a/interface/HZZ4LRooPdfs.h +++ b/interface/HZZ4LRooPdfs.h @@ -10,6 +10,8 @@ #include "TDirectory.h" #include "TH2F.h" +#include // for ROOT_VERSION + #include #include #include @@ -136,6 +138,10 @@ class RooqqZZPdf_v2 : public RooAbsPdf { RooqqZZPdf_v2(const RooqqZZPdf_v2& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new RooqqZZPdf_v2(*this,newname); } inline ~RooqqZZPdf_v2() override { } + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override; +#endif protected: @@ -299,6 +305,10 @@ class RooggZZPdf_v2 : public RooAbsPdf { RooggZZPdf_v2(const RooggZZPdf_v2& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new RooggZZPdf_v2(*this,newname); } inline ~RooggZZPdf_v2() override { } + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + std::unique_ptr compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const override; +#endif protected: diff --git a/interface/ProcessNormalization.h b/interface/ProcessNormalization.h index 24916cc3792..a3444c5d2e7 100644 --- a/interface/ProcessNormalization.h +++ b/interface/ProcessNormalization.h @@ -4,6 +4,8 @@ #include #include "RooListProxy.h" +#include "CombineCodegenImpl.h" + //_________________________________________________ /* BEGIN_HTML @@ -19,7 +21,6 @@ class ProcessNormalization : public RooAbsReal { ProcessNormalization(const char *name, const char *title, double nominal=1) ; ProcessNormalization(const char *name, const char *title, RooAbsReal &nominal) ; ProcessNormalization(const ProcessNormalization &other, const char *newname = 0) ; - ~ProcessNormalization() override ; TObject * clone(const char *newname) const override { return new ProcessNormalization(*this, newname); } @@ -29,6 +30,15 @@ class ProcessNormalization : public RooAbsReal { void addOtherFactor(RooAbsReal &factor) ; void dump() const ; + COMBINE_DECLARE_TRANSLATE; + + double nominalValue() const { return nominalValue_; } + std::vector const &logKappa() const { return logKappa_; } + RooArgList const &thetaList() const { return thetaList_; } + std::vector> const &logAsymmKappa() const { return logAsymmKappa_; } + RooArgList const &asymmThetaList() const { return asymmThetaList_; } + RooArgList const &otherFactorList() const { return otherFactorList_; } + protected: Double_t evaluate() const override; @@ -51,4 +61,6 @@ class ProcessNormalization : public RooAbsReal { ClassDefOverride(ProcessNormalization,1) // Process normalization interpolator }; +COMBINE_DECLARE_CODEGEN_IMPL(ProcessNormalization); + #endif diff --git a/interface/VerticalInterpHistPdf.h b/interface/VerticalInterpHistPdf.h index f1172b6f6d0..ccfc34aad26 100644 --- a/interface/VerticalInterpHistPdf.h +++ b/interface/VerticalInterpHistPdf.h @@ -10,7 +10,8 @@ #include "TH1.h" #include "SimpleCacheSentry.h" #include "FastTemplate_Old.h" -#include + +#include "CombineCodegenImpl.h" class FastVerticalInterpHistPdf; class FastVerticalInterpHistPdf2Base; @@ -34,7 +35,12 @@ class VerticalInterpHistPdf : public RooAbsPdf { const RooArgList& coefList() const { return _coefList ; } const RooRealProxy &x() const { return _x; } + + double smoothRegion() const { return _smoothRegion; } + int smoothAlgo() const { return _smoothAlgo; } + friend class FastVerticalInterpHistPdf; + protected: RooRealProxy _x; RooListProxy _funcList ; // List of component FUNCs @@ -84,6 +90,11 @@ class FastVerticalInterpHistPdfBase : public RooAbsPdf { /// Must be public, for serialization struct Morph { FastTemplate sum; FastTemplate diff; }; + std::vector const &morphs() const { return _morphs; } + + double smoothRegion() const { return _smoothRegion; } + int smoothAlgo() const { return _smoothAlgo; } + friend class FastVerticalInterpHistPdf2Base; protected: RooRealProxy _x; @@ -127,21 +138,22 @@ class FastVerticalInterpHistPdf : public FastVerticalInterpHistPdfBase { FastVerticalInterpHistPdf() : FastVerticalInterpHistPdfBase() {} FastVerticalInterpHistPdf(const char *name, const char *title, const RooRealVar &x, const RooArgList& funcList, const RooArgList& coefList, Double_t smoothRegion=1., Int_t smoothAlgo=1) : FastVerticalInterpHistPdfBase(name,title,RooArgSet(x),funcList,coefList,smoothRegion,smoothAlgo), - _x("x","Independent variable",this,const_cast(x)), - _cache(), _cacheNominal(), _cacheNominalLog() {} + _x("x","Independent variable",this,const_cast(x)) {} FastVerticalInterpHistPdf(const FastVerticalInterpHistPdf& other, const char* name=0) : FastVerticalInterpHistPdfBase(other, name), _x("x",this,other._x), _cache(other._cache), _cacheNominal(other._cacheNominal), _cacheNominalLog(other._cacheNominalLog) {} TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf(*this,newname) ; } - ~FastVerticalInterpHistPdf() override {} Double_t evaluate() const override ; Bool_t hasCache() const { return _cache.size() > 0; } Bool_t isCacheReady() const { return _cache.size() > 0 && _init; } + FastHisto const &cacheNominal() const { return _cacheNominal; } + friend class FastVerticalInterpHistPdfV; friend class FastVerticalInterpHistPdf2; + protected: RooRealProxy _x; @@ -186,14 +198,12 @@ class FastVerticalInterpHistPdf2D : public FastVerticalInterpHistPdfBase { FastVerticalInterpHistPdfBase(name,title,RooArgSet(x,y),funcList,coefList,smoothRegion,smoothAlgo), _x("x","Independent variable",this,const_cast(x)), _y("y","Independent variable",this,const_cast(y)), - _conditional(conditional), - _cache(), _cacheNominal(), _cacheNominalLog() {} + _conditional(conditional) {} FastVerticalInterpHistPdf2D(const FastVerticalInterpHistPdf2D& other, const char* name=0) : FastVerticalInterpHistPdfBase(other, name), _x("x",this,other._x), _y("y",this,other._y), _conditional(other._conditional), _cache(other._cache), _cacheNominal(other._cacheNominal), _cacheNominalLog(other._cacheNominalLog) {} TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf2D(*this,newname) ; } - ~FastVerticalInterpHistPdf2D() override {} const RooRealVar & x() const { return dynamic_cast(_x.arg()); } const RooRealVar & y() const { return dynamic_cast(_y.arg()); } @@ -204,6 +214,8 @@ class FastVerticalInterpHistPdf2D : public FastVerticalInterpHistPdfBase { Bool_t hasCache() const { return _cache.size() > 0; } Bool_t isCacheReady() const { return _cache.size() > 0 && _init; } + FastHisto2D const &cacheNominal() const { return _cacheNominal; } + friend class FastVerticalInterpHistPdf2D2; protected: RooRealProxy _x, _y; @@ -245,8 +257,15 @@ class FastVerticalInterpHistPdf2Base : public RooAbsPdf { virtual void setActiveBins(unsigned int bins) {} bool cacheIsGood() const { return _sentry.good() && _initBase; } + + double smoothRegion() const { return _smoothRegion; } + int smoothAlgo() const { return _smoothAlgo; } + /// Must be public, for serialization typedef FastVerticalInterpHistPdfBase::Morph Morph; + + std::vector const &morphs() const { return _morphs; } + protected: RooListProxy _coefList ; // List of coefficients Double_t _smoothRegion; @@ -302,7 +321,6 @@ class FastVerticalInterpHistPdf2 : public FastVerticalInterpHistPdf2Base { _cache(other._cache), _cacheNominal(other._cacheNominal), _cacheNominalLog(other._cacheNominalLog) {} explicit FastVerticalInterpHistPdf2(const FastVerticalInterpHistPdf& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf2(*this,newname) ; } - ~FastVerticalInterpHistPdf2() override {} const RooRealVar & x() const { return dynamic_cast(_x.arg()); } @@ -311,6 +329,10 @@ class FastVerticalInterpHistPdf2 : public FastVerticalInterpHistPdf2Base { FastHisto const& cache() const { return _cache; } + COMBINE_DECLARE_TRANSLATE; + + FastHisto const &cacheNominal() const { return _cacheNominal; } + friend class FastVerticalInterpHistPdf2V; protected: RooRealProxy _x; @@ -358,7 +380,6 @@ class FastVerticalInterpHistPdf2D2 : public FastVerticalInterpHistPdf2Base { _cache(other._cache), _cacheNominal(other._cacheNominal), _cacheNominalLog(other._cacheNominalLog) {} explicit FastVerticalInterpHistPdf2D2(const FastVerticalInterpHistPdf2D& other, const char* name=0) ; TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf2D2(*this,newname) ; } - ~FastVerticalInterpHistPdf2D2() override {} const RooRealVar & x() const { return dynamic_cast(_x.arg()); } const RooRealVar & y() const { return dynamic_cast(_y.arg()); } @@ -369,6 +390,10 @@ class FastVerticalInterpHistPdf2D2 : public FastVerticalInterpHistPdf2Base { Double_t evaluate() const override ; + COMBINE_DECLARE_TRANSLATE; + + FastHisto2D const &cacheNominal() const { return _cacheNominal; } + protected: RooRealProxy _x, _y; bool _conditional; @@ -400,14 +425,12 @@ class FastVerticalInterpHistPdf3D : public FastVerticalInterpHistPdfBase { _x("x","Independent variable",this,const_cast(x)), _y("y","Independent variable",this,const_cast(y)), _z("z","Independent variable",this,const_cast(z)), - _conditional(conditional), - _cache(), _cacheNominal(), _cacheNominalLog() {} + _conditional(conditional) {} FastVerticalInterpHistPdf3D(const FastVerticalInterpHistPdf3D& other, const char* name=0) : FastVerticalInterpHistPdfBase(other, name), _x("x",this,other._x), _y("y",this,other._y), _z("z",this,other._z), _conditional(other._conditional), _cache(other._cache), _cacheNominal(other._cacheNominal), _cacheNominalLog(other._cacheNominalLog) {} TObject* clone(const char* newname) const override { return new FastVerticalInterpHistPdf3D(*this,newname) ; } - ~FastVerticalInterpHistPdf3D() override {} const RooRealVar & x() const { return dynamic_cast(_x.arg()); } const RooRealVar & y() const { return dynamic_cast(_y.arg()); } @@ -441,5 +464,7 @@ class FastVerticalInterpHistPdf3D : public FastVerticalInterpHistPdfBase { }; +COMBINE_DECLARE_CODEGEN_IMPL(FastVerticalInterpHistPdf2); +COMBINE_DECLARE_CODEGEN_IMPL(FastVerticalInterpHistPdf2D2); #endif diff --git a/scripts/fitRooFitAD.py b/scripts/fitRooFitAD.py new file mode 100644 index 00000000000..314e21a1d59 --- /dev/null +++ b/scripts/fitRooFitAD.py @@ -0,0 +1,52 @@ +import numpy as np +import ROOT +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("--input", "-i", help="input ws file") +parser.add_argument("--backend", help='set evaluation backend, default is "combine"') +parser.add_argument("--write-debug-macro", action="store_true", help='write debug macro in case of the "codegen" backend') + + +args = parser.parse_args() + +with ROOT.TFile.Open(args.input) as f: + ws = f.Get("w") + +global_observables = ws.set("globalObservables") +constrain = ws.set("nuisances") + +pdf = ws["model_s"] +data = ws["data_obs"] + +# To use the evaluation backends of RooFit +ROOT.gInterpreter.Declare("#include ") + +ROOT.Combine.nllBackend_ = args.backend + +# Change verbosity +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Minimization) +ROOT.RooMsgService.instance().getStream(1).removeTopic(ROOT.RooFit.Fitting) + +ROOT.gInterpreter.Declare("#include ") + +nll = pdf.createNLL(data, Constrain=constrain, GlobalObservables=global_observables, Offset="initial") + +ROOT.gInterpreter.Declare( + """ +void writeDebugMacro(RooAbsReal & real, std::string const &name) +{ + static_cast(real).writeDebugMacro(name); +} +""" +) + +if args.backend == "codegen" and args.write_debug_macro: + ROOT.writeDebugMacro(nll, "codegen") + +cfg = ROOT.RooMinimizer.Config() +minim = ROOT.RooMinimizer(nll, cfg) +minim.setEps(1.0) +minim.setStrategy(0) +minim.minimize("Minuit2", "") +minim.save().Print() diff --git a/src/AsymPow.cc b/src/AsymPow.cc index fb7ba951c13..1a3ed897e20 100644 --- a/src/AsymPow.cc +++ b/src/AsymPow.cc @@ -20,13 +20,6 @@ AsymPow::AsymPow(const AsymPow &other, const char *newname) : theta_("theta",this,other.theta_) { } -AsymPow::~AsymPow() {} - -TObject *AsymPow::clone(const char *newname) const -{ - return new AsymPow(*this,newname); -} - Double_t AsymPow::evaluate() const { return RooFit::Detail::MathFuncs::asymPow(theta_, kappaLow_, kappaHigh_); } diff --git a/src/CombineCodegenImpl.cxx b/src/CombineCodegenImpl.cxx new file mode 100644 index 00000000000..405a9695003 --- /dev/null +++ b/src/CombineCodegenImpl.cxx @@ -0,0 +1,213 @@ +#include "../interface/CombineCodegenImpl.h" + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) + +#include "../interface/AsymPow.h" +#include "../interface/ProcessNormalization.h" +#include "../interface/VerticalInterpHistPdf.h" + +#include + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +namespace RooFit { +namespace Experimental { +# define CODEGEN_IMPL(CLASS_NAME) void codegenImpl(CLASS_NAME &arg0, CodegenContext &ctx) +# define ARG_VAR auto &arg = arg0; +#else +# define CODEGEN_IMPL(CLASS_NAME) void CLASS_NAME::translate(RooFit::Detail::CodeSquashContext &ctx) const +# define ARG_VAR auto &arg = *this; +#endif + + +CODEGEN_IMPL(AsymPow) { + ARG_VAR; + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::asymPow", arg.theta(), arg.kappaLow(), arg.kappaHigh())); +} + +CODEGEN_IMPL(ProcessNormalization) { + ARG_VAR; + + std::vector logAsymmKappaLow; + std::vector logAsymmKappaHigh; + logAsymmKappaLow.reserve(arg.logAsymmKappa().size()); + logAsymmKappaHigh.reserve(arg.logAsymmKappa().size()); + for (auto [lo, hi] : arg.logAsymmKappa()) { + logAsymmKappaLow.push_back(lo); + logAsymmKappaHigh.push_back(hi); + } + + ctx.addResult(&arg, + ctx.buildCall("RooFit::Detail::MathFuncs::processNormalization", + arg.nominalValue(), + arg.thetaList().size(), + arg.asymmThetaList().size(), + arg.otherFactorList().size(), + arg.thetaList(), + arg.logKappa(), + arg.asymmThetaList(), + logAsymmKappaLow, + logAsymmKappaHigh, + arg.otherFactorList())); +} + +CODEGEN_IMPL(FastVerticalInterpHistPdf2) { + ARG_VAR; + + if (arg.smoothAlgo() < 0) { + throw std::runtime_error("We only support smoothAlgo >= 0"); + } + + RooRealVar const &xVar = arg.x(); + + int numBins = xVar.numBins(); + + std::vector nominalVec(numBins); + std::vector widthVec(numBins); + std::vector morphsVecSum; + std::vector morphsVecDiff; + + auto const &cacheNominal = arg.cacheNominal(); + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + widthVec[i] = cacheNominal.GetWidth(i); + } + + std::size_t nCoefs = arg.coefList().size(); + + morphsVecSum.reserve(numBins * nCoefs); + morphsVecDiff.reserve(numBins * nCoefs); + auto const &morphs = arg.morphs(); + for (unsigned int j = 0; j < nCoefs; ++j) { + for (int i = 0; i < numBins; ++i) { + morphsVecSum.push_back(morphs[j].sum[i]); + morphsVecDiff.push_back(morphs[j].diff[i]); + } + } + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + } + + // The bin index part + // We also have to assert that x is uniformely binned! + if (!dynamic_cast(&xVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + double xLow = xVar.getMin(); + double xHigh = xVar.getMax(); + std::string binIdx = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, xVar, numBins); + + std::string arrName = ctx.getTmpVarName(); + + std::stringstream code; + code << "double " << arrName << "[" << numBins << "];\n"; + code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2", + numBins, + nCoefs, + arg.coefList(), + nominalVec, + widthVec, + morphsVecSum, + morphsVecDiff, + arg.smoothRegion(), + arrName) + + ";\n"; + + ctx.addToCodeBody(code.str(), true); + ctx.addResult(&arg, arrName + "[" + binIdx + "]"); +} + +CODEGEN_IMPL(FastVerticalInterpHistPdf2D2) { + ARG_VAR; + + if (arg.smoothAlgo() < 0) { + throw std::runtime_error("We only support smoothAlgo >= 0"); + } + + if (!arg.conditional()) { + throw std::runtime_error("We only support conditional == true"); + } + + RooRealVar const &xVar = arg.x(); + RooRealVar const &yVar = arg.y(); + + // We also have to assert that x and y are uniformely binned! + if (!dynamic_cast(&xVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + if (!dynamic_cast(&yVar.getBinning())) { + throw std::runtime_error("We only support uniform binning!"); + } + + auto const &cacheNominal = arg.cacheNominal(); + + int numBinsX = cacheNominal.binX(); + int numBinsY = cacheNominal.binY(); + int numBins = numBinsY * numBinsY; + + std::vector nominalVec(numBins); + std::vector widthVec(numBins); + std::vector morphsVecSum; + std::vector morphsVecDiff; + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + widthVec[i] = cacheNominal.GetWidth(i); + } + + std::size_t nCoefs = arg.coefList().size(); + + morphsVecSum.reserve(numBins * nCoefs); + morphsVecDiff.reserve(numBins * nCoefs); + auto const &morphs = arg.morphs(); + for (unsigned int j = 0; j < nCoefs; ++j) { + for (int i = 0; i < numBins; ++i) { + morphsVecSum.push_back(morphs[j].sum[i]); + morphsVecDiff.push_back(morphs[j].diff[i]); + } + } + + for (int i = 0; i < numBins; ++i) { + nominalVec[i] = cacheNominal.GetBinContent(i); + } + + // The bin index part + double xLow = xVar.getMin(); + double xHigh = xVar.getMax(); + std::string binIdxX = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", xLow, xHigh, arg.x(), numBinsX); + double yLow = yVar.getMin(); + double yHigh = yVar.getMax(); + std::string binIdxY = ctx.buildCall("RooFit::Detail::MathFuncs::getUniformBinning", yLow, yHigh, arg.y(), numBinsY); + + std::stringstream binIdx; + binIdx << "(" << binIdxY << " + " << yVar.numBins() << " * " << binIdxX << ")"; + + std::string arrName = ctx.getTmpVarName(); + + std::stringstream code; + code << "double " << arrName << "[" << (numBinsX * numBinsY) << "];\n"; + code << ctx.buildCall("RooFit::Detail::MathFuncs::fastVerticalInterpHistPdf2D2", + numBinsX, + numBinsY, + nCoefs, + arg.coefList(), + nominalVec, + widthVec, + morphsVecSum, + morphsVecDiff, + arg.smoothRegion(), + arrName) + + ";\n"; + + ctx.addToCodeBody(code.str(), true); + ctx.addResult(&arg, arrName + "[" + binIdx.str() + "]"); +} + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,35,0) +} // namespace RooFit +} // namespace Experimental +#endif + +#endif // ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) diff --git a/src/HZZ4LRooPdfs.cc b/src/HZZ4LRooPdfs.cc index 613031da83b..302638513b3 100644 --- a/src/HZZ4LRooPdfs.cc +++ b/src/HZZ4LRooPdfs.cc @@ -4,7 +4,7 @@ #include "RooAbsReal.h" #include "RooRealVar.h" #include "RooVoigtian.h" -//#include "RooComplex.h" +#include "RooGenericPdf.h" #include "RooMath.h" #include "RooAbsCategory.h" #include @@ -2126,6 +2126,15 @@ Double_t RooqqZZPdf_v2::evaluate() const } +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +std::unique_ptr RooqqZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const +{ + return RooGenericPdf{GetName(), + "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3])))+(.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))+(.5+.5*TMath::Erf((x[0]-x[11])/x[12]))*(x[14]/(1+exp((x[0]-x[11])/x[13])) )", + {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9, *a10, *a11, *a12, *a13} + }.compileForNormSet(normSet, ctx); +} +#endif /************RooggZZPdf_v2***********/ @@ -2189,6 +2198,16 @@ Double_t RooggZZPdf_v2::evaluate() const return ZZ; } +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,32,0) +std::unique_ptr RooggZZPdf_v2::compileForNormSet(RooArgSet const &normSet, RooFit::Detail::CompileContext & ctx) const +{ + return RooGenericPdf{GetName(), + "(.5+.5*TMath::Erf((x[0]-x[1])/x[2]))*(x[4]/(1+exp((x[0]-x[1])/x[3]))) + (.5+.5*TMath::Erf((x[0]-x[5])/x[6]))*(x[8]/(1+exp((x[0]-x[5])/x[7]))+x[10]/(1+exp((x[0]-x[5])/x[9])))", + {*m4l, *a0, *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8, *a9} + }.compileForNormSet(normSet, ctx); +} +#endif + //// ------- RooVBFZZPdf ------ diff --git a/src/ProcessNormalization.cc b/src/ProcessNormalization.cc index 5baa229b3a1..b1fbb936388 100644 --- a/src/ProcessNormalization.cc +++ b/src/ProcessNormalization.cc @@ -15,14 +15,10 @@ ProcessNormalization::ProcessNormalization(const char *name, const char *title, { } -ProcessNormalization::ProcessNormalization(const char *name, const char *title, RooAbsReal &nominal) : - RooAbsReal(name,title), - nominalValue_(1.0), - thetaList_("thetaList", "List of nuisances for symmetric kappas", this), - asymmThetaList_("asymmThetaList", "List of nuisances for asymmetric kappas", this), - otherFactorList_("otherFactorList", "Other multiplicative terms", this) -{ - otherFactorList_.add(nominal); +ProcessNormalization::ProcessNormalization(const char *name, const char *title, RooAbsReal &nominal) + : ProcessNormalization{name, title, 1.0} +{ + otherFactorList_.add(nominal); } ProcessNormalization::ProcessNormalization(const ProcessNormalization &other, const char *newname) : @@ -36,8 +32,6 @@ ProcessNormalization::ProcessNormalization(const ProcessNormalization &other, co { } -ProcessNormalization::~ProcessNormalization() {} - void ProcessNormalization::addLogNormal(double kappa, RooAbsReal &theta) { if (kappa != 0.0 && kappa != 1.0) { logKappa_.push_back(std::log(kappa)); diff --git a/src/VerticalInterpHistPdf.cc b/src/VerticalInterpHistPdf.cc index ff9ef2ec65a..df2ac8ee5ed 100644 --- a/src/VerticalInterpHistPdf.cc +++ b/src/VerticalInterpHistPdf.cc @@ -10,6 +10,8 @@ #include "RooMsgService.h" #include "RooAbsData.h" +#include "RooUniformBinning.h" + //#define TRACE_CALLS #ifdef TRACE_CALLS #include "../interface/ProfilingTools.h"