include/root/RooNDKeysPdf.h

0001 /*****************************************************************************
0002  * Project: RooFit                                                           *
0003  * Package: RooFitModels                                                     *
0004  *    File: $Id: RooNDKeysPdf.h 44368 2012-05-30 15:38:44Z axel $
0005  * Authors:                                                                  *
0006  *   Max Baak, CERN, mbaak@cern.ch *
0007  *                                                                           *
0008  * Copyright (c) 2000-2005, Regents of the University of California          *
0009  *                          and Stanford University. All rights reserved.    *
0010  *                                                                           *
0011  * Redistribution and use in source and binary forms,                        *
0012  * with or without modification, are permitted according to the terms        *
0013  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
0014  *****************************************************************************/
0015 #ifndef ROO_NDKEYS_PDF
0016 #define ROO_NDKEYS_PDF
0017
0018 #include "RooAbsPdf.h"
0019 #include "RooRealProxy.h"
0020 #include "RooSetProxy.h"
0021 #include "RooDataSet.h"
0022 #include "RooListProxy.h"
0023 #include "TH1.h"
0024 #include "TAxis.h"
0025 #include "TVectorD.h"
0026 #include "TMatrixD.h"
0027 #include "TMatrixDSym.h"
0028 #include <map>
0029 #include <vector>
0030 #include <string>
0031
0032 class RooRealVar;
0033 class RooArgList;
0034 class RooArgSet;
0035 class RooChangeTracker;
0036
0037 class VecVecDouble : public std::vector<std::vector<double> >  { } ;
0038 class VecTVecDouble : public std::vector<TVectorD> { } ;
0039 typedef std::pair<Int_t, VecVecDouble::iterator > iiPair;
0040 typedef std::vector< iiPair > iiVec;
0041 typedef std::pair<Int_t, VecTVecDouble::iterator > itPair;
0042 typedef std::vector< itPair > itVec;
0043
0044 class RooNDKeysPdf : public RooAbsPdf {
0045
0046 public:
0047
0048   enum Mirror {NoMirror, MirrorLeft, MirrorRight, MirrorBoth,
0049                MirrorAsymLeft, MirrorAsymLeftRight,
0050                MirrorAsymRight, MirrorLeftAsymRight,
0051                MirrorAsymBoth };
0052
0053   RooNDKeysPdf() = default;
0054
0055   RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
0056                TString options = "ma", double rho = 1, double nSigma = 3, bool rotate = true,
0057                bool sortInput = true);
0058
0059   RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist, TString options = "ma",
0060                double rho = 1, double nSigma = 3, bool rotate = true, bool sortInput = true);
0061
0062   RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
0063                const TVectorD &rho, TString options = "ma", double nSigma = 3, bool rotate = true,
0064                bool sortInput = true);
0065
0066   RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const RooDataSet &data,
0067                const RooArgList &rhoList, TString options = "ma", double nSigma = 3, bool rotate = true,
0068                bool sortInput = true);
0069
0070   RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, const TH1 &hist,
0071                const RooArgList &rhoList, TString options = "ma", double nSigma = 3, bool rotate = true,
0072                bool sortInput = true);
0073
0074   RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, const RooDataSet &data, Mirror mirror = NoMirror,
0075                double rho = 1, double nSigma = 3, bool rotate = true, bool sortInput = true);
0076
0077   RooNDKeysPdf(const char *name, const char *title, RooAbsReal &x, RooAbsReal &y, const RooDataSet &data,
0078                TString options = "ma", double rho = 1.0, double nSigma = 3, bool rotate = true,
0079                bool sortInput = true);
0080
0081   RooNDKeysPdf(const RooNDKeysPdf& other, const char* name=nullptr);
0082   ~RooNDKeysPdf() override;
0083
0084   TObject* clone(const char* newname) const override { return new RooNDKeysPdf(*this,newname); }
0085
0086   Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=nullptr) const override ;
0087   double analyticalIntegral(Int_t code, const char* rangeName=nullptr) const override ;
0088
0089   TMatrixD getWeights(const int &k) const;
0090
0091   struct BoxInfo {
0092     bool filled;
0093     bool netFluxZ;
0094     double nEventsBW;
0095     double nEventsBMSW;
0096     std::vector<double> xVarLo, xVarHi;
0097     std::vector<double> xVarLoM3s, xVarLoP3s, xVarHiM3s, xVarHiP3s;
0098     std::map<Int_t,bool> bpsIdcs;
0099     std::vector<Int_t> sIdcs;
0100     std::vector<Int_t> bIdcs;
0101     std::vector<Int_t> bmsIdcs;
0102   } ;
0103
0104 protected:
0105
0106   RooListProxy _varList ;
0107   RooListProxy _rhoList;
0108
0109   double evaluate() const override;
0110
0111   void createPdf(bool firstCall, RooDataSet const& data);
0112   void setOptions();
0113   void initialize(RooDataSet const& data);
0114   void loadDataSet(bool firstCall, RooDataSet const& data);
0115   void mirrorDataSet();
0116   void loadWeightSet(RooDataSet const& data);
0117   void calculateShell(BoxInfo *bi) const;
0118   void calculatePreNorm(BoxInfo *bi) const;
0119   void sortDataIndices(BoxInfo *bi = nullptr);
0120   void calculateBandWidth();
0121   double gauss(std::vector<double> &x, std::vector<std::vector<double>> &weights) const;
0122   void loopRange(std::vector<double> &x, std::vector<Int_t> &indices) const;
0123   void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const;
0124   RooDataSet *createDatasetFromHist(const RooArgList &varList, const TH1 &hist) const;
0125   void updateRho() const;
0126   void checkInitWeights() const {
0127     if (_weights == &_weights0 || _weights == &_weights1)
0128       return;
0129     const_cast<RooNDKeysPdf*>(this)->calculateBandWidth();
0130   }
0131
0132   mutable TString _options;
0133   double _widthFactor;
0134   double _nSigma;
0135
0136   bool _fixedShape{false};
0137   bool _mirror{false};
0138   bool _debug{false};   //!
0139   bool _verbose{false}; //!
0140
0141   Int_t _nDim{0};
0142   Int_t _nEvents{0};
0143   Int_t _nEventsM{0};
0144   double _nEventsW{0.};
0145   double _d{0.};
0146   double _n{0.};
0147
0148   // cached info on variable
0149   std::vector<std::vector<double> > _dataPts;
0150   std::vector<TVectorD> _dataPtsR;
0151   std::vector<std::vector<double> > _weights0;            // Plain weights
0152   std::vector<std::vector<double> > _weights1;            // Weights for adaptive kernels
0153   std::vector<std::vector<double> >* _weights{nullptr};   //! Weights to be used. Points either to _weights0 or _weights1
0154
0155   std::vector<itVec> _sortTVIdcs; //!
0156
0157   mutable std::vector<double> _rho;
0158   RooArgSet _dataVars;
0159   mutable std::vector<double> _x; // Cache for x values
0160   std::vector<double> _x0, _x1, _x2;
0161   std::vector<double> _mean, _sigma;
0162   std::vector<double> _xDatLo, _xDatHi;
0163   std::vector<double> _xDatLo3s, _xDatHi3s;
0164
0165   bool _netFluxZ{false};
0166   double _nEventsBW{0.};
0167   double _nEventsBMSW{0.};
0168   std::vector<double> _xVarLo, _xVarHi;
0169   std::vector<double> _xVarLoM3s, _xVarLoP3s, _xVarHiM3s, _xVarHiP3s;
0170   std::map<Int_t,bool> _bpsIdcs;
0171   std::map<Int_t, bool> _ibNoSort;
0172   std::vector<Int_t> _sIdcs;
0173   std::vector<Int_t> _bIdcs;
0174   std::vector<Int_t> _bmsIdcs;
0175
0176   // Data for analytical integrals:
0177   mutable std::map<std::pair<std::string,int>,BoxInfo*> _rangeBoxInfo ;
0178   mutable BoxInfo _fullBoxInfo ;
0179
0180   std::vector<Int_t> _idx;
0181   double _minWeight{0.};
0182   double _maxWeight{0.};
0183   std::map<Int_t,double> _wMap;
0184
0185   TMatrixDSym* _covMat{nullptr};
0186   TMatrixDSym* _corrMat{nullptr};
0187   TMatrixD* _rotMat{nullptr};
0188   TVectorD* _sigmaR{nullptr};
0189   TVectorD* _dx{nullptr};
0190   double _sigmaAvgR{0.};
0191
0192   bool _rotate;
0193   bool _sortInput;
0194   Int_t _nAdpt;
0195
0196   RooChangeTracker *_tracker{nullptr}; //
0197
0198   ClassDefOverride(RooNDKeysPdf, 3) // General N-dimensional non-parametric kernel estimation p.d.f
0199 };
0200
0201 #endif