From 0e4929c6358d65e6c56d322691158a2088243edb Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Tue, 23 Jul 2024 15:15:58 +0200 Subject: [PATCH] modify stats --- .aux/test_with_lcg | 2 +- ostap/core/core.py | 4 +- ostap/fitting/dataset.py | 48 +- ostap/frames/frames.py | 33 +- ostap/stats/statvars.py | 188 +++++-- ostap/trees/cuts.py | 43 +- ostap/trees/trees.py | 8 + source/CMake_OSTAP.cmake | 3 +- source/include/Ostap/Covariance.h | 251 ++++++++++ source/include/Ostap/Moments.h | 9 +- source/include/Ostap/StatEntity.h | 11 +- source/include/Ostap/StatVar.h | 564 +++++---------------- source/include/Ostap/WStatEntity.h | 15 +- source/src/Covariance.cpp | 278 +++++++++++ source/src/StatEntity.cpp | 11 + source/src/StatVar.cpp | 753 ++++++----------------------- source/src/WStatEntity.cpp | 10 + source/src/dict/Ostap.hh | 1 + 18 files changed, 1084 insertions(+), 1148 deletions(-) create mode 100644 source/include/Ostap/Covariance.h create mode 100644 source/src/Covariance.cpp diff --git a/.aux/test_with_lcg b/.aux/test_with_lcg index e739d831..2d559b64 100755 --- a/.aux/test_with_lcg +++ b/.aux/test_with_lcg @@ -4,5 +4,5 @@ CMTCONFIG=$2 source /cvmfs/sft.cern.ch/lcg/views/${LCG}/${CMTCONFIG}/setup.sh source build/INSTALL/thisostap.sh cd build -ctest -N && cmake .. -DCMAKE_INSTALL_PREFIX=./INSTALL/ && ctest -j8 -R '(stat|frame)' --output-on-failure +ctest -N && cmake .. -DCMAKE_INSTALL_PREFIX=./INSTALL/ && ctest -j8 --output-on-failure diff --git a/ostap/core/core.py b/ostap/core/core.py index d278b95f..60c21ffb 100644 --- a/ostap/core/core.py +++ b/ostap/core/core.py @@ -1058,7 +1058,7 @@ def split_string_respect ( text , separators = ',;:' , strip = True ) : ## strip items if required if strip : items = [ item.strip() for item in items ] ## remove empty items - return [ item for item in items if item ] + return tuple ( item for item in items if item ) # ============================================================================= ## split string using separators: @@ -1088,7 +1088,7 @@ def split_string ( line , ## strip items if required if strip : items = [ i.strip() for i in items ] ## remove empty items - return [ item for item in items if item ] + return tuple ( item for item in items if item ) diff --git a/ostap/fitting/dataset.py b/ostap/fitting/dataset.py index 82637967..f39bace6 100644 --- a/ostap/fitting/dataset.py +++ b/ostap/fitting/dataset.py @@ -37,7 +37,7 @@ list_types , sequence_types ) from ostap.math.base import islong from ostap.fitting.variables import valid_formula, make_formula -import ostap.trees.cuts +from ostap.trees.cuts import expression_types, vars_and_cuts import ostap.fitting.roocollections import ostap.fitting.printable import ROOT, random, math, sys, ctypes @@ -1143,6 +1143,14 @@ def ds_var_minmax ( dataset , var , cuts = '' , delta = 0.0 ) : >>> mn,mx = data.vminmax('pt') >>> mn,mx = data.vminmax('pt','y>3') """ + + assert isinstance ( var , expression_types ) , 'Invalid expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid expression!' + + var = str ( var ) + cuts = str ( cuts ).strip() + + if isinstance ( var , ROOT.RooAbsReal ) : var = var.GetName() if cuts : s = dataset.statVar ( var , cuts ) else : s = dataset.statVar ( var ) @@ -1188,8 +1196,12 @@ def _ds_has_entry_ ( dataset , selection , *args ) : >>> dataset.hasEntry ( 'pt>100' , 0, 1000 ) ## >>> dataset.hasEntry ( 'pt>100' , 'fit_range' ) ## >>> dataset.hasEntry ( 'pt>100' , 'fit_range' , 0 , 1000 ) ## - - see Ostap.StatVar.hasEntru + - see Ostap.StatVar.hasEntry """ + + assert isinstance ( selection , expression_types ) , 'Invalid expression!' + selection = str ( selection ).strip() + result = Ostap.StatVar.hasEntry ( dataset , selection , *args ) return True if result else False @@ -1853,24 +1865,12 @@ def _ds_table_0_ ( dataset , if not valid_pointer ( varset ) : logger.error('Invalid dataset') return '' - - if isinstance ( variables , str ) : - variables = split_string ( variables , var_separators , strip = True , respect_groups = True ) - - if 1 == len ( variables ) : variables = variables [0] - - if isinstance ( variables , str ) : - - if variables in varset : - vars = [ variables ] - else : - vars = list ( dataset.branches ( variables ) ) - - elif variables : vars = [ i.GetName() for i in varset if i.GetName() in variables ] - else : vars = [ i.GetName() for i in varset ] + + if variables : vars , cuts = vars_and_cuts ( variables , cuts ) + else : vars , cuts = [ v.name for v in varset ] , str ( cuts ).strip() + + vars = [ i.GetName() for i in varset if i.GetName() in vars ] - # - _vars = [] stat = dataset.statVars ( vars , cuts , first , last ) @@ -2067,11 +2067,12 @@ def _ds_table_1_ ( dataset , """Print data set as table """ - if isinstance ( variables , str ) : - variables = split_string ( variables , var_separators , strip = True , respect_groups = True ) - + + if variables : vars , cuts = vars_and_cuts ( variables , cuts ) + else : vars , cuts = [ v.name for v in varset ] , str ( cuts ).strip() + _vars = [] - vvars = tuple ( sorted ( variables ) ) + vvars = tuple ( sorted ( vars ) ) stat = dataset.statVars ( vvars , cuts , first , last ) for v in stat : s = stat [ v ] @@ -2087,7 +2088,6 @@ def _ds_table_1_ ( dataset , _vars.append ( r ) _vars.sort() - tt = dataset.GetTitle() if not title : diff --git a/ostap/frames/frames.py b/ostap/frames/frames.py index 52e51756..399fadc2 100644 --- a/ostap/frames/frames.py +++ b/ostap/frames/frames.py @@ -345,37 +345,6 @@ def frame_statCov ( frame , ## possible types of expressions expression_types = string_types + ( ROOT.TCut , ) # ============================================================================== -## Prepare the arguments: variable names and cuts -# @code -# vars_lst, cuts, input_string = vars_and_cuts ( 'x', 'y<0' ) -# vars_lst, cuts, input_string = vars_and_cuts ( ('x','y','z') , 'y<0' ) -# @endcode -def vars_and_cuts ( expressions , cuts ) : - """Prepare the arguments: variable names and cuts - >>> vars_lst, cuts, input_string = vars_and_cuts ( 'x', 'y<0' ) - >>> vars_lst, cuts, input_string = vars_and_cuts ( ('x','y','z') , 'y<0' ) - """ - - ## single string as input - input_string = False - if isinstance ( expressions , expression_types ) : - input_string = True - expressions = [ expressions ] - - assert expressions and all ( s and isinstance ( s , expression_types ) for s in expressions ) , \ - "Invalid expression(s) : %s" % str ( expressions ) - - exprs = tuple ( str(s).strip() for s in expressions ) - assert exprs and all ( exprs ) , "Invalid expression(s): %s" % str ( exprs ) - - assert isinstance ( cuts , expression_types ) or not cuts , \ - 'Invaild type of cuts: %s' % str ( cuts ) - - cuts = str ( cuts ) if cuts else '' - - return exprs, cuts, input_string - -# ================================================================================== ## helper functionm that defines expression and cuts # @code # frame = ... @@ -387,7 +356,7 @@ def _fr_helper_ ( frame , expressions , cuts = '' ) : >>> current , vexpr , cexpr = _fr_helper_ ( frame , 'x*x' , 'z<0' ) """ - exprs, cuts, input_string = vars_and_cuts ( expressions , cuts ) + exprs, cuts, input_string = SV.vars_and_cuts ( expressions , cuts ) ## Frame/Tree ? if isinstance ( frame , ROOT.TTree ) : frame = DataFrame ( frame ) diff --git a/ostap/stats/statvars.py b/ostap/stats/statvars.py index eebfbf82..51502054 100644 --- a/ostap/stats/statvars.py +++ b/ostap/stats/statvars.py @@ -72,13 +72,17 @@ 'data_power_mean' , ## get the (weighted) power mean 'data_lehmer_mean' , ## get the (weighted) Lehmer mean 'data_decorate' , ## technical function to decorate the class - ) + 'expression_types' , ## valid types for expressions/cuts/weights +) # ============================================================================= from builtins import range -from ostap.core.core import Ostap, rootException -from ostap.core.ostap_types import string_types, num_types +from ostap.math.base import isequal, iszero +from ostap.core.core import Ostap, rootException, strings, WSE +from ostap.core.ostap_types import string_types, integer_types, num_types +from ostap.trees.cuts import expression_types, vars_and_cuts import ostap.stats.moment -import ostap.logger.table as T +import ostap.logger.table as T +import ROOT # ============================================================================= # logging # ============================================================================= @@ -91,7 +95,6 @@ ## @var QEXACT # use it as threshold for exact/slow vs approximate/fast quantile calcualtion QEXACT = 10000 - # ============================================================================= ## get the moment of order 'order' relative to 'center' # @code @@ -107,7 +110,15 @@ def data_get_moment ( data , order , center , expression , cuts = '' , *args ) : >>> print data.get_moment ( 3 , 0.0 , 'mass' , 'pt>1' ) ## ditto - see Ostap::StatVar::get_moment """ - assert isinstance ( order , int ) and 0<= order , 'Invalid order %s' % order + assert isinstance ( order , integer_types ) and 0<= order , 'Invalid order %s' % order + assert isinstance ( center , num_types ) , 'Invalid center!' + + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + with rootException() : return StatVar.get_moment ( data , order , @@ -130,7 +141,13 @@ def data_moment ( data , order , expression , cuts = '' , *args ) : >>> print data.moment ( 3 , 'mass' , 'pt>1' ) ## ditto - see Ostap::StatVar::moment """ - assert isinstance ( order , int ) and 0<= order , 'Invalid order %s' % order + assert isinstance ( order , integer_types ) and 0 <= order , 'Invalid order %s' % order + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + with rootException() : return StatVar.moment ( data , order , @@ -158,7 +175,13 @@ def data_central_moment ( data , order , expression , cuts = '' , *args ) : - Uncertainty is calculated with O(1/n^2) precision - see Ostap::StatVar::central_moment """ - assert isinstance ( order , int ) and 0<= order , 'Invalid order %s' % order + assert isinstance ( order , integer_types ) and 0 <= order , 'Invalid order %s' % order + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + with rootException() : return StatVar.central_moment ( data , order , @@ -166,7 +189,7 @@ def data_central_moment ( data , order , expression , cuts = '' , *args ) : cuts , *args ) # ============================================================================== -## Get the statsitocs from data +## Get the statistics from data # @code # statobj = Ostap.Math.MinValue() # data = ... @@ -186,18 +209,96 @@ def data_get_stat ( data , statobj , expression , cuts = '' , *args ) : assert isinstance ( statobj , ( Ostap.Math.Statistic , Ostap.Math.WStatistic ) ) ,\ 'get_object: invalid statobj type: %s' % type ( statobj ) - import ROOT - if ( not cuts ) and isinstance ( data , ROOT.TTree ) and isinstance ( statobj , Ostap.Math.Statistic ) : + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + + if isinstance ( data , ROOT.TTree ) and isinstance ( statobj , Ostap.Math.Statistic ) and not cuts : with rootException() : sc = StatVar.the_moment ( data , statobj , expression , *args ) assert sc.isSuccess() , 'Error %s from StatVar::the_moment' % sc return statobj - + with rootException() : sc = StatVar.the_moment ( data , statobj , expression , cuts , *args ) assert sc.isSuccess() , 'Error %s from StatVar::the_moment' % sc return statobj +# ============================================================================== +## Get the statistics from data +# @code +# data = ... +# result = data_statistics ( data , 'x+y' , 'pt>1' ) +# @encode +# @see Ostap::StatEntity +# @see Ostap::WStatEntity +# @see Ostap::StatVar::statVar +def data_statistics ( data , expressions , cuts = '' , *args ) : + """Get statistics from data + >>> data = ... + >>> result = data_statistics ( data , 'x/y+z' , '0>> data = ... + >>> result = data_covariance ( data , 'x/y+z' , 'qq' , 'u>0') + - see Ostap.Math.Covariance + - see Ostap.Math.WCovariance + - see Ostap.StatVar.statCov + """ + assert isinstance ( expression1 , expression_types ) , 'Invalid type of expression1!' + assert isinstance ( expression2 , expression_types ) , 'Invalid type of expression2!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression1 = str ( expression1 ) + expression2 = str ( expression2 ) + cuts = str ( cuts ).strip() + + with rootException() : + if isinstance ( data , ROOT.TTree ) and not cuts : + return StatVar.statCov ( data , expression1 , expression2 , *args ) + else : + return StatVar.statCov ( data , expression1 , expression2 , cuts , *args ) + # ============================================================================= ## Get harmonic mean over the data # @code @@ -218,14 +319,11 @@ def data_harmonic_mean ( data , expression , cuts = '' , *args ) : - see `Ostap.Math::HarmonicMean` - see `Ostap.Math::WHarmonicMean` - see `Ostap.statVar.the_moment` - """ - - import ROOT - - if ( not cuts ) and isinstance ( data , ROOT.TTree ) : + """ + if isinstance ( data , ROOT.TTree ) and not cuts : statobj = Ostap.Math.HarmonicMean() return data_get_stat ( data , statobj , expression , '' , *args ) - + statobj = Ostap.Math.WHarmonicMean() return data_get_stat ( data , statobj , expression , cuts , *args ) @@ -248,13 +346,10 @@ def data_geometric_mean ( data , expression , cuts = '' , *args ) : - see `Ostap.Math::HarmonicMean` - see `Ostap.statVar.the_moment` """ - - import ROOT - if ( not cuts ) and isinstance ( data , ROOT.TTree ) : statobj = Ostap.Math.GeometricMean() return data_get_stat ( data , statobj , expression , '' , *args ) - + statobj = Ostap.Math.WGeometricMean() return data_get_stat ( data , statobj , expression , cuts , *args ) @@ -279,17 +374,13 @@ def data_power_mean ( data , p , expression , cuts = '' , *args ) : - see `Ostap.Math::WPowerMean` - see `Ostap.statVar.the_moment` """ - import ROOT - from ostap.core.ostap_types import num_types - from ostap.math.base import isequal, iszero - assert isinstance ( p , num_types ) , 'Invalid p-parameter type: %s' % type ( p ) if p == -1 or isequal ( p , -1. ) : return data_harmonic_mean ( data , expression , cuts , *args ) elif p == 0 or iszero ( p ) : return data_geometric_mean ( data , expression , cuts , *args ) elif p == 1 or isequal ( p , 1. ) : return data_arithmetic_mean ( data , expression , cuts , *args ) - if ( not cuts ) and isinstance ( data , ROOT.TTree ) : + if isinstance ( data , ROOT.TTree ) and not cuts : statobj = Ostap.Math.PowerMean ( p ) return data_get_stat ( data , statobj , expression , '' , *args ) @@ -317,10 +408,8 @@ def data_arithmetic_mean ( data , expression , cuts = '' , *args ) : - see `Ostap.Math::WArithmeticMean` - see `Ostap.statVar.the_moment` """ - import ROOT - from ostap.math.base import isequal, iszero - if ( not cuts ) and isinstance ( data , ROOT.TTree ) : + if isinstance ( data , ROOT.TTree ) and not cuts : statobj = Ostap.Math.ArithmeticMean ( ) return data_get_stat ( data , statobj , expression , '' , *args ) @@ -348,9 +437,6 @@ def data_lehmer_mean ( data , p , expression , cuts = '' , *args ) : - see `Ostap.Math::WLehmerMean` - see `Ostap.statVar.the_moment` """ - import ROOT - from ostap.core.ostap_types import num_types - from ostap.math.base import isequal, iszero assert isinstance ( p , num_types ) , 'Invalid p-parameter!' @@ -364,7 +450,6 @@ def data_lehmer_mean ( data , p , expression , cuts = '' , *args ) : statobj = Ostap.Math.WLehmerMean ( p ) return data_get_stat ( data , statobj , expression , cuts , *args ) - # ============================================================================= ## Get the moments or order order as Ostap::Math::(W)Moment_ # @code @@ -382,7 +467,13 @@ def data_the_moment ( data , order , expression , cuts = '' , *args ) : """ assert isinstance ( order , int ) and 0 <= order , 'Invalid order %s' % order - if ( not cuts ) and isinstance ( data , ROOT.TTree ) : + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + + if isinstance ( data , ROOT.TTree ) and not cuts : M = Ostap.Math. Moment_(order) moment = M () with rootException() : @@ -474,6 +565,13 @@ def data_skewness ( data , expression , cuts = '' , *args ) : >>> print data.skewness ( 'mass' , 'pt>1' ) ## ditto - see Ostap::StatVar::skewness """ + + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + with rootException() : return StatVar.skewness ( data , expression , cuts , *args ) @@ -514,6 +612,12 @@ def data_quantile ( data , q , expression , cuts = '' , exact = QEXACT , *args """ assert isinstance ( q , float ) and 0 < q < 1 , 'Invalid quantile:%s' % q + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + if exact is True : ## exact slow algorithm with rootException() : @@ -555,6 +659,12 @@ def data_interval ( data , qmin , qmax , expression , cuts = '' , exact = QEXAC assert isinstance ( qmin , float ) and 0 < qmin < 1 , 'Invalid quantile-1:%s' % qmin assert isinstance ( qmax , float ) and 0 < qmax < 1 , 'Invalid quantile-2:%s' % qmax + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + qmin, qmax = min ( qmin , qmax ) , max ( qmin, qmax ) if exact is True : @@ -621,6 +731,12 @@ def data_quantiles ( data , quantiles , expression , cuts = '' , exact = QEXACT >>> print data.quantiles ( 10 , 'mass' , 'pt>1' ) ## deciles! - see Ostap::StatVar::quantile """ + assert isinstance ( expression , expression_types ) , 'Invalid type of expression!' + assert isinstance ( cuts , expression_types ) , 'Invalid type of cuts/weight!' + + expression = str ( expression ) + cuts = str ( cuts ).strip() + if isinstance ( quantiles , float ) and 0 < quantiles < 1 : quantiles = [ quantiles ] elif isinstance ( quantiles , int ) and 1 < quantiles : @@ -923,8 +1039,6 @@ def _qi_str_ ( o ) : return "QInterval([%.5g,%.5g],n=%s)" % ( o.interval.low , Interval .__repr__ = _i_str_ QInterval .__str__ = _qi_str_ QInterval .__repr__ = _qi_str_ - - # ============================================================================= diff --git a/ostap/trees/cuts.py b/ostap/trees/cuts.py index c94626c1..2be7fac7 100755 --- a/ostap/trees/cuts.py +++ b/ostap/trees/cuts.py @@ -13,10 +13,13 @@ __version__ = "$Revision$" __author__ = "Vanya BELYAEV Ivan.Belyaev@itep.ru" __date__ = "2011-06-07" -__all__ = () ## nothing to import +__all__ = ( + 'expression_types' , ## valid types for expressions/selections/weights + 'vars_and_cuts' , ## helper routibe to treat expressions +) # ============================================================================= -from ostap.core.core import cpp, VE, hID, dsID from ostap.core.ostap_types import num_types, string_types +from ostap.core.core import cpp, VE, hID, dsID, split_string from ostap.utils.utils import balanced import ROOT # ============================================================================= @@ -27,6 +30,42 @@ else : logger = getLogger( __name__ ) # ============================================================================= logger.debug( 'Some useful decorations for ROOT.TCut objects') +# ============================================================================= +## types for expressions and cuts +expression_types = string_types + ( ROOT.TCut , ) +# ============================================================================= +## Prepare the arguments: variable names and cuts +# @code +# vars_lst, cuts, input_string = vars_and_cuts ( 'x', 'y<0' ) +# vars_lst, cuts, input_string = vars_and_cuts ( ('x','y','z') , 'y<0' ) +# @endcode +def vars_and_cuts ( expressions , cuts ) : + """Prepare the arguments: variable names and cuts + >>> vars_lst, cuts, input_string = vars_and_cuts ( 'x', 'y<0' ) + >>> vars_lst, cuts, input_string = vars_and_cuts ( ('x','y','z') , 'y<0' ) + """ + + ## single string as input + input_string = False + if isinstance ( expressions , expression_types ) : + expressions = split_string ( str ( expressions ) , strip = True , respect_groups = True ) + input_string = 1 == len ( expressions ) + + assert expressions and all ( s and isinstance ( s , expression_types ) for s in expressions ) , \ + "Invalid expression(s) : %s" % str ( expressions ) + + exprs = tuple ( str(s).strip() for s in expressions ) + exprs = tuple ( s for s in expressions if s ) + assert exprs and all ( exprs ) , "Invalid expression(s): %s" % str ( exprs ) + + assert isinstance ( cuts , expression_types ) or not cuts , \ + 'Invaild type of cuts: %s' % str ( cuts ) + + cuts = str ( cuts ).strip() if cuts else '' + + return exprs , cuts, input_string + + # ============================================================================= # ROOT.TCut # ============================================================================= diff --git a/ostap/trees/trees.py b/ostap/trees/trees.py index d990bdf9..e27aaba8 100755 --- a/ostap/trees/trees.py +++ b/ostap/trees/trees.py @@ -57,6 +57,8 @@ # ============================================================================= + + # ============================================================================= ## check validity/emptiness of TTree/TChain # require non-zero poniter and non-empty Tree/Chain @@ -827,6 +829,9 @@ def _stat_cov_ ( tree , Use only subset of events >>> stat1 , stat2 , cov2 , len = tree.statCov( 'x' , 'y' , 'z>0' , 100 , 10000 ) """ + + + """ import ostap.math.linalg stat1 = Ostap.WStatEntity () stat2 = Ostap.WStatEntity () @@ -850,7 +855,10 @@ def _stat_cov_ ( tree , *args ) return stat1 , stat2 , cov2, length + """ + + ROOT.TTree . statCov = _stat_cov_ ROOT.TChain . statCov = _stat_cov_ diff --git a/source/CMake_OSTAP.cmake b/source/CMake_OSTAP.cmake index 97f0dd82..74b67d23 100644 --- a/source/CMake_OSTAP.cmake +++ b/source/CMake_OSTAP.cmake @@ -13,8 +13,9 @@ add_library(ostap SHARED src/format.cpp src/BreitWigner.cpp src/ChebyshevApproximation.cpp src/Choose.cpp - src/Combine.cpp src/Chi2Fit.cpp + src/Combine.cpp + src/Covariance.cpp src/Dalitz.cpp src/DalitzIntegrator.cpp src/DataFrameActions.cpp diff --git a/source/include/Ostap/Covariance.h b/source/include/Ostap/Covariance.h new file mode 100644 index 00000000..5b7fe0b1 --- /dev/null +++ b/source/include/Ostap/Covariance.h @@ -0,0 +1,251 @@ +// ============================================================================ +#ifndef OSTAP_COVARIANCE_H +#define OSTAP_COVARIANCE_H 1 +// ============================================================================ +// Include files +// ============================================================================ +// STD & STL +// ============================================================================ +#include +// ============================================================================ +// Ostap +// ============================================================================ +#include "Ostap/StatEntity.h" +#include "Ostap/WStatEntity.h" +#include "Ostap/SymmetricMatrixTypes.h" +// ============================================================================ +namespace Ostap +{ + // ========================================================================== + namespace Math + { + // ======================================================================== + /** @class Covariance + * Counter for two variables which also counts the covariance + * @see Ostap::StatEntity + */ + class Covariance + { + // ====================================================================== + public: + // ====================================================================== + /// the actual type of counter + typedef Ostap::StatEntity Counter ; + typedef Ostap::SymMatrix2x2 Matrix ; + // ====================================================================== + public: + // ====================================================================== + /// default constructor + Covariance () = default ; + // ====================================================================== + // constructor from two counters and the correlation coefficients + Covariance + ( const Counter& c1 , + const Counter& c2 , + const double corr = 0 ) ; + // ====================================================================== + public: + // ====================================================================== + /// the first counter + const Counter& counter1 () const { return m_cnt1 ; } + /// the second counter + const Counter& counter2 () const { return m_cnt2 ; } + // ====================================================================== + /// get the moment \f$ \sum_i (x_i - \bar{x} ) ( y_i - \bar{y}) \f$ + inline double cov2m () const { return m_cov2m ; } + /// get the true covarinace + inline double covariance () const { return empty() ? 0.0 : m_cov2m / n() ; } + /// get the correlation coefficient + double correlation () const ; + // ====================================================================== + /// number of entries + inline unsigned long long n () const { return m_cnt1.n () ; } + /// effective number of entries + inline unsigned long long nEff () const { return m_cnt1.nEff () ; } + /// empty ? + inline bool empty () const { return m_cnt1.empty () ; } + // ====================================================================== + public: + // ====================================================================== + /// add two values to the counters + inline Covariance& add + ( const double x , + const double y ) ; + /// add another counter + Covariance& add ( const Covariance& right ) ; + // ====================================================================== + /// add another counter + inline Covariance& operator+= ( const Covariance& right ) + { return add ( right ) ; } + // ====================================================================== + public: + // ====================================================================== + /// add x,y + Covariance& update + ( const double x , + const double y ) { return add ( x , y ) ; } + /// add another counter + Covariance& update ( const Covariance& right ) { return add ( right ) ; } + // ====================================================================== + /// reset counters + void reset () ; + // ====================================================================== + private: + // ====================================================================== + Ostap::StatEntity m_cnt1 { } ; + Ostap::StatEntity m_cnt2 { } ; + double m_cov2m { 0 } ; + // ====================================================================== + }; + // ======================================================================== + /// add two values to the counters + inline Covariance& + Covariance::add + ( const double x , + const double y ) + { + const unsigned long long nn = n () ; + if ( nn ) + { + const double xA = m_cnt1.mean() ; + const double yA = m_cnt2.mean() ; + m_cov2m += ( x - xA ) * nn * ( y - yA ) / ( nn + 1 ) ; + } + // + m_cnt1 += x ; + m_cnt2 += y ; + // + return *this ; + } + // ======================================================================== + /// external operator for addition of two covariance objects + inline Covariance operator+ ( Covariance a , const Covariance& b ) + { a += b ; return a ; } + // ======================================================================== + /// get the covariance matrix + Covariance::Matrix covariance ( const Covariance& ) ; + // ======================================================================== + /// get the correlation matrix + Covariance::Matrix correlation ( const Covariance& ) ; + // ======================================================================== + /** @class WCovariance + * Counter for two variables which also counts the covariance + * @see Ostap::WStatEntity + */ + class WCovariance + { + // ====================================================================== + public: + // ====================================================================== + /// the actual type of counter + typedef Ostap::WStatEntity Counter ; + typedef Ostap::SymMatrix2x2 Matrix ; + // ====================================================================== + public: + // ====================================================================== + /// default constructor + WCovariance () = default ; + // ====================================================================== + // constructor from two counters and the correlation coefficients + WCovariance + ( const Counter& c1 , + const Counter& c2 , + const double corr = 0 ) ; + // ====================================================================== + public: + // ====================================================================== + /// the first counter + const Counter& counter1 () const { return m_cnt1 ; } + /// the second counter + const Counter& counter2 () const { return m_cnt2 ; } + // ====================================================================== + /// get the moment \f$ \sum_i (x_i - \bar{x} ) ( y_i - \bar{y}) \f$ + inline double cov2m () const { return m_cov2m ; } + /// get the true covarinace + inline double covariance () const { return empty() ? 0.0 : m_cov2m / w () ; } + /// get the correlation coefficient + double correlation () const ; + // ====================================================================== + /// number of entries + inline unsigned long long n () const { return m_cnt1.n () ; } + /// effective number of entries + inline double nEff () const { return m_cnt1.nEff () ; } + /// empty ? + inline bool empty () const { return m_cnt1.empty () ; } + /// sum of weights + inline double w () const { return m_cnt1.sumw () ; } + /// sum of weights + inline double sumw () const { return m_cnt1.sumw () ; } + // ====================================================================== + public: + // ====================================================================== + /// add two values to the counters + inline WCovariance& add + ( const double x , + const double y , + const double w = 1 ) ; + /// add another counter + WCovariance& add ( const WCovariance& right ) ; + // ====================================================================== + /// add another counter + inline WCovariance& operator+= ( const WCovariance& right ) + { return add ( right ) ; } + // ====================================================================== + public: + // ====================================================================== + /// add x,y + WCovariance& update + ( const double x , + const double y , + const double w = 1 ) { return add ( x , y , w ) ; } + /// add another counter + WCovariance& update ( const WCovariance& right ) { return add ( right ) ; } + // ====================================================================== + /// reset counters + void reset () ; + // ====================================================================== + private: + // ====================================================================== + Ostap::WStatEntity m_cnt1 { } ; + Ostap::WStatEntity m_cnt2 { } ; + double m_cov2m { 0 } ; + // ====================================================================== + }; + // ======================================================================== + /// add two values to the counters + inline WCovariance& + WCovariance::add + ( const double x , + const double y , + const double w ) + { + const double ww = sumw () ; + if ( ww && w ) + { + const double xA = m_cnt1.mean() ; + const double yA = m_cnt2.mean() ; + m_cov2m += ( x - xA ) * ww * w * ( y - yA ) / ( ww + w ) ; + } + // + m_cnt1.add ( x , w ) ; + m_cnt2.add ( y , w ) ; + // + return *this ; + } + // ======================================================================== + /// external operator for addition of two covariance objects + inline WCovariance operator+ ( WCovariance a , const WCovariance& b ) + { a += b ; return a ; } + // ======================================================================== + /// get the covariance matrix + Covariance::Matrix covariance ( const WCovariance& ) ; + // ======================================================================== + /// get the correlation matrix + Covariance::Matrix correlation ( const WCovariance& ) ; + // ======================================================================== + } // The end of nameapace Ostap::Math + // ========================================================================== +} // The end namespace Ostap +// ============================================================================ +#endif // OSTAP_COVARIANCE_H +// ============================================================================ diff --git a/source/include/Ostap/Moments.h b/source/include/Ostap/Moments.h index 61c39ab0..401a2783 100644 --- a/source/include/Ostap/Moments.h +++ b/source/include/Ostap/Moments.h @@ -2095,12 +2095,7 @@ namespace Ostap Moment_<0> m_cnt {} ; // ====================================================================== } ; - - - - - - + // ======================================================================== /** @class Moments * Collection of static functions dealing with moments @@ -2561,7 +2556,7 @@ namespace Ostap // ====================================================================== } ; // The end of class Ostap::Math::Moments - // ======================================================================== + // ======================================================================== } // The end of namespace Ostap::Math // ========================================================================== } // The end of namespace Ostap diff --git a/source/include/Ostap/StatEntity.h b/source/include/Ostap/StatEntity.h index 21120702..d434858a 100644 --- a/source/include/Ostap/StatEntity.h +++ b/source/include/Ostap/StatEntity.h @@ -70,6 +70,8 @@ namespace Ostap // ====================================================================== /// number of entries unsigned long long n () const { return m_n ; } + /// effective number of entries + unsigned long long nEff () const { return m_n ; } /// mean value double mu () const { return m_mu ; } /// the second central moment/dispersion/variance @@ -81,6 +83,8 @@ namespace Ostap // ====================================================================== public: // derived quantities & aliases // ====================================================================== + /// empty ? + bool empty () const { return 0 == m_n ; } /// number of entries unsigned long long nEntries () const { return m_n ; } /// variance @@ -313,6 +317,8 @@ namespace Ostap // ======================================================================== /// reset the counters void reset () ; + /// swap two counters + void swap ( StatEntity& right ) ; /// representation as string std::string toString () const; /// printout to std::ostream @@ -348,7 +354,10 @@ namespace Ostap // ========================================================================== /// conversion to string inline std::string to_string ( const StatEntity& s ) { return s.toString() ; } - // ========================================================================== + // ========================================================================== + /// swap two counters + inline void swap ( StatEntity& a , StatEntity& b ) { a.swap( b ) ; } + // ========================================================================= } // end of namespace Ostap // ============================================================================ // The END diff --git a/source/include/Ostap/StatVar.h b/source/include/Ostap/StatVar.h index 431471ba..98f94fd3 100644 --- a/source/include/Ostap/StatVar.h +++ b/source/include/Ostap/StatVar.h @@ -12,13 +12,14 @@ // ============================================================================= class TTree ; // ROOT class TChain ; // ROOT -class TCut ; // ROOT class RooAbsData ; // RooFit // ============================================================================= // Ostap // ============================================================================ +#include "Ostap/StatEntity.h" #include "Ostap/WStatEntity.h" #include "Ostap/ValueWithError.h" +#include "Ostap/Covariance.h" #include "Ostap/SymmetricMatrixTypes.h" #include "Ostap/DataFrame.h" #include "Ostap/StatusCode.h" @@ -27,8 +28,8 @@ class RooAbsData ; // RooFit // ============================================================================ template class TMatrixTSym ; // ROOT // ============================================================================ -namespace Ostap { namespace Math { class Statistic ; } } -namespace Ostap { namespace Math { class WStatistic ; } } +namespace Ostap { namespace Math { class Statistic ; } } +namespace Ostap { namespace Math { class WStatistic ; } } // ============================================================================ namespace Ostap { @@ -36,7 +37,6 @@ namespace Ostap /** @class StatVar Ostap/StatVar.h * Helper class to get statistical * infomation about the variable/expression - * * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2013-10-13 */ @@ -45,11 +45,15 @@ namespace Ostap public: // ======================================================================== /// the actual type for statistic - typedef Ostap::WStatEntity Statistic ; + typedef Ostap::WStatEntity Statistic ; /// the actual type for vector of statistic - typedef std::vector Statistics ; + typedef std::vector Statistics ; /// variable names - typedef std::vector Names ; + typedef std::vector Names ; + /// covariance + typedef Ostap::Math::Covariance Covariance ; + /// (weighted) covariance + typedef Ostap::Math::WCovariance WCovariance ; // ======================================================================== public: // ======================================================================== @@ -58,7 +62,9 @@ namespace Ostap */ struct Interval { - Interval ( const double l = 0 , const double h = 0 ) + Interval + ( const double l = 0 , + const double h = 0 ) : low ( std::min ( l , h ) ) , high ( std::max ( l , h ) ) {} @@ -75,7 +81,9 @@ namespace Ostap */ struct Quantile { - Quantile ( const double q = 0 , const unsigned long n = 0 ) + Quantile + ( const double q = 0 , + const unsigned long n = 0 ) : quantile ( q ) , nevents ( n ) {} @@ -92,8 +100,9 @@ namespace Ostap */ struct Quantiles { - Quantiles ( const std::vector& q = std::vector() , - const unsigned long n = 0 ) + Quantiles + ( const std::vector& q = std::vector() , + const unsigned long n = 0 ) : quantiles ( q ) , nevents ( n ) {} @@ -110,8 +119,9 @@ namespace Ostap */ struct QInterval { - QInterval ( const Interval& i = Interval() , - const unsigned long n = 0 ) + QInterval + ( const Interval& i = Interval () , + const unsigned long n = 0 ) : interval ( i ) , nevents ( n ) {} @@ -131,26 +141,6 @@ namespace Ostap // ======================================================================== public: // ======================================================================== - /** build statistic for the expression - * @param tree (INPUT) the tree - * @param expression (INPUT) the expression - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * - * @code - * tree = ... - * stat = tree.statVar ( 'S_sw' ) - * @endcode - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2013-10-13 - */ - static Statistic statVar - ( TTree* tree , - const std::string& expression , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** build statistic for the expression * @param tree (INPUT) the tree * @param expression (INPUT) the expression @@ -160,16 +150,15 @@ namespace Ostap * * @code * tree = ... - * stat = tree.statVar( 'S_sw' ,'pt>1000') + * stat = tree.statVar( 'S_sw' ) * @endcode * * @author Vanya BELYAEV Ivan.Belyaev@itep.ru * @date 2013-10-13 */ - static Statistic statVar + static Ostap::StatEntity statVar ( TTree* tree , const std::string& expression , - const std::string& cuts , const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== @@ -191,7 +180,7 @@ namespace Ostap static Statistic statVar ( TTree* tree , const std::string& expression , - const TCut& cuts , + const std::string& cuts , const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== @@ -210,19 +199,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** check if there is at least one entry that satisfies criteria - * @param tree (INPUT) the tree - * @param cuts (INPUT) criteria - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ - static bool hasEntry - ( TTree* tree , - const TCut& cuts , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** check if there is at least one entry that satisfies criteria * @param data (INPUT) data * @param cuts (INPUT) criteria @@ -236,19 +212,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** check if there is at least one entry that satisfies criteria - * @param data (INPUT) data - * @param cuts (INPUT) criteria - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ - static bool hasEntry - ( const RooAbsData* data , - const TCut& cuts , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** check if there is at least one entry that satisfies criteria * @param data (INPUT) data * @param cuts (INPUT) criteria @@ -264,21 +227,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** check if there is at least one entry that satisfies criteria - * @param data (INPUT) data - * @param cuts (INPUT) criteria - * @param cut_range (INPUT) cut range - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ - static bool hasEntry - ( const RooAbsData* data , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== public: // ======================================================================== /** build statistic for the expressions @@ -317,25 +265,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** build statistic for the expressions - * @param tree (INPUT) the tree - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - * @param cuts (INPUT) the selection criteria - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-11-04 - */ - static unsigned long statVars - ( TTree* tree , - std::vector& result , - const Names& expressions , - const TCut& cuts , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== public: // ======================================================================== /** build statistic for the expression @@ -382,28 +311,6 @@ namespace Ostap const unsigned long last = LAST ) { return statVar( data , expression , cuts , "" , first , last ) ; } // ======================================================================== - /** build statistic for the expression - * @param data (INPUT) the data - * @param expression (INPUT) the expression - * @param cuts (INPUT) the selection - * @param first (INPUT) the frist event to process - * @param last (INPUT) process till "last"-event - * - * @code - * data = ... - * stat = data.statVar( 'S_sw' , 'pt>10') - * @endcode - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2015-02-15 - */ - static Statistic statVar - ( const RooAbsData* data , - const std::string& expression , - const TCut& cuts , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** build statistic for the expression * @param data (INPUT) the data * @param expression (INPUT) the expression @@ -429,31 +336,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** build statistic for the expression - * @param data (INPUT) the data - * @param expression (INPUT) the expression - * @param cuts (INPUT) the selection - * @param cutrange (INPUT) the cut-range - * @param first (INPUT) the frist event to process - * @param last (INPUT) process till "last"-event - * - * @code - * data = ... - * cut = TCut ( ... ) - * stat = data.statVar( 'S_sw' , cut ) - * @endcode - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2015-02-15 - */ - static Statistic statVar - ( const RooAbsData* data , - const std::string& expression , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== public: // ======================================================================== /** build statistic for the expressions @@ -494,25 +376,6 @@ namespace Ostap const unsigned long last = LAST ) { return statVars ( data , result , expressions , cuts , std::string() , first , last ) ; } // ======================================================================== - /** build statistic for the expressions - * @param data (INPUT) input data - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - * @param cuts (INPUT) the selection - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2021-06-04 - */ - static unsigned long statVars - ( const RooAbsData* data , - std::vector& result , - const Names& expressions , - const TCut& cuts , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** build statistic for the expressions * @param data (INPUT) input data * @param result (UPDATE) the output statistics for specified expressions @@ -534,26 +397,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** build statistic for the expressions - * @param data (INPUT) input data - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - # @param cuts (INPUT) the selection - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2021-06-04 - */ - static unsigned long statVars - ( const RooAbsData* data , - std::vector& result , - const Names& expressions , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== public: // ======================================================================== /** build statistic for the expression @@ -573,78 +416,52 @@ namespace Ostap ( FrameNode frame , const std::string& expression , const std::string& cuts = "" ) ; - // ======================================================================== + // ======================================================================== public: // ======================================================================== - /** calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 + /** get the number of equivalent entries + * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ + * @param tree (INPUT) the tree + * @param cuts (INPUT) selection criteria + * @param first (INPUT) the first event to process + * @param last (INPUT) the last event to process + * @return number of equivalent entries */ - static unsigned long statCov - ( TTree* tree , - const std::string& exp1 , - const std::string& exp2 , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , + static double nEff + ( TTree& tree , + const std::string& cuts = "" , const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection criteria - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 + /** get the number of equivalent entries + * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ + * @param tree (INPUT) the tree + * @param cuts (INPUT) selection crietria + * @param first (INPUT) the first event to process + * @param last (INPUT) the last event to process + * @param cut_range (INPUT) cut range + * @return number of equivalent entries */ - static unsigned long statCov - ( TTree* tree , - const std::string& exp1 , - const std::string& exp2 , - const std::string& cuts , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; + static double nEff + ( const RooAbsData& tree , + const std::string& cuts = "" , + const std::string& cut_range = "" , + const unsigned long first = 0 , + const unsigned long last = LAST ) ; // ======================================================================== - /** calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection criteria - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 + /** get the number of equivalent entries + * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ + * @param tree (INPUT) the frame + * @param cuts (INPUT) selection crietria/weight + * @param first (INPUT) the first event to process + * @param last (INPUT) the last event to process + * @return number of equivalent entries */ - static unsigned long statCov - ( TTree* tree , - const std::string& exp1 , - const std::string& exp2 , - const TCut& cuts , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; + static double nEff + ( FrameNode frame , + const std::string& cuts = "" ) ; + // ======================================================================== + public: // ======================================================================== /** calculate the covariance of several expressions * @param tree (INPUT) the inpout tree @@ -665,24 +482,41 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** calculate the covariance of several expressions + public: // new methods + // ======================================================================== + /** calculate the covariance of two expressions * @param tree (INPUT) the inpout tree - * @param vars (INPUT) expressions - * @param cuts (INPUT) the selection criteria - * @param stats (UPDATE) the statistics - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2023-02-28 + * @param exp1 (INPUT) the first expression + * @param exp2 (INPUT) the second expression + * @return Covariance + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2024-07-22- */ - static unsigned long statCov - ( TTree* tree , - const std::vector& vars , - const TCut& cuts , - std::vector& stats , - TMatrixTSym& cov2 , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; + static Covariance statCov + ( TTree* tree , + const std::string& exp1 , + const std::string& exp2 , + const unsigned long first = 0 , + const unsigned long last = LAST ) ; + // ======================================================================== + /** calculate the covariance of two expressions + * @param tree (INPUT) the inpout tree + * @param exp1 (INPUT) the first expression + * @param exp2 (INPUT) the second expression + * @param cuts (INPUT) cuts/weigths expression + * @return Covariance + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2024-07-22 + */ + static WCovariance statCov + ( TTree* tree , + const std::string& exp1 , + const std::string& exp2 , + const std::string& cuts , + const unsigned long first = 0 , + const unsigned long last = LAST ) ; + // ======================================================================== + public: // ======================================================================== /** calculate the covariance of several expressions * @param tree (INPUT) the inpout tree @@ -703,54 +537,6 @@ namespace Ostap // ======================================================================== public: // ======================================================================== - /** calculate the covariance of two expressions - * @param tree (INPUT) the input tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 - */ - static unsigned long statCov - ( const RooAbsData* tree , - const std::string& exp1 , - const std::string& exp2 , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const std::string& cut_range = "" , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== - /** calculate the covariance of two expressions - * @param tree (INPUT) the input tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 - */ - static unsigned long statCov - ( const RooAbsData* tree , - const std::string& exp1 , - const std::string& exp2 , - const std::string& cuts , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const std::string& cut_range = "" , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** calculate the covariance of several expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions @@ -772,27 +558,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** calculate the covariance of several expressions - * @param tree (INPUT) the inpout tree - * @param vars (INPUT) expressions - * @param cuts (INPUT) the selection criteria - * @param stats (UPDATE) the statistics - * @param cov2 (UPDATE) the covariance matrix - * @param cut_range (INPUT) range - * @return number of processed events - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2023-02-28 - */ - static unsigned long statCov - ( const RooAbsData* tree , - const std::vector& vars , - const TCut& cuts , - std::vector& stats , - TMatrixTSym& cov2 , - const std::string& cut_range = "" , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** calculate the covariance of several expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions @@ -815,90 +580,53 @@ namespace Ostap public: // ======================================================================== /** calculate the covariance of two expressions - * @param frame (INPUT) data frame + * @param tree (INPUT) the input tree * @param exp1 (INPUT) the first expresiion * @param exp2 (INPUT) the second expresiion - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix + * @param cuts (INPUT) selection/weight expression * @return number of processed events * * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-06-18 + * @date 2014-03-27 */ - static unsigned long statCov - ( FrameNode data , + static WCovariance statCov + ( const RooAbsData* tree , const std::string& exp1 , const std::string& exp2 , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 ) ; + const std::string& cuts , + const std::string& cut_range = "" , + const unsigned long first = 0 , + const unsigned long last = LAST ) ; + // ======================================================================== + public: // ======================================================================== /** calculate the covariance of two expressions * @param frame (INPUT) data frame * @param exp1 (INPUT) the first expresiion * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection/weight - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * + * @return Covariancce * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-06-18 + * @date 2024-07-22 */ - static unsigned long statCov - ( FrameNode data , - const std::string& exp1 , - const std::string& exp2 , - const std::string& cuts , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 ) ; - // ======================================================================== - public: + static Covariance statCov + ( FrameNode data , + const std::string& exp1 , + const std::string& exp2 ) ; // ======================================================================== - /** get the number of equivalent entries - * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ - * @param tree (INPUT) the tree - * @param cuts (INPUT) selection criteria - * @param first (INPUT) the first event to process - * @param last (INPUT) the last event to process - * @return number of equivalent entries - */ - static double nEff - ( TTree& tree , - const std::string& cuts = "" , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== - /** get the number of equivalent entries - * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ - * @param tree (INPUT) the tree - * @param cuts (INPUT) selection crietria - * @param first (INPUT) the first event to process - * @param last (INPUT) the last event to process - * @param cut_range (INPUT) cut range - * @return number of equivalent entries - */ - static double nEff - ( const RooAbsData& tree , - const std::string& cuts = "" , - const std::string& cut_range = "" , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== - /** get the number of equivalent entries - * \f$ n_{eff} \equiv = \frac{ (\sum w)^2}{ \sum w^2} \f$ - * @param tree (INPUT) the frame - * @param cuts (INPUT) selection crietria/weight - * @param first (INPUT) the first event to process - * @param last (INPUT) the last event to process - * @return number of equivalent entries + /** calculate the covariance of two expressions + * @param frame (INPUT) data frame + * @param exp1 (INPUT) the first expresiion + * @param exp2 (INPUT) the second expresiion + * @param cuts (INPUT) selection/weight expresiion + * @return Covariancce + * @author Vanya BELYAEV Ivan.Belyaev@itep.ru + * @date 2024-07-22 */ - static double nEff - ( FrameNode frame , - const std::string& cuts = "" ) ; + static WCovariance statCov + ( FrameNode data , + const std::string& exp1 , + const std::string& exp2 , + const std::string& cuts ) ; // ======================================================================== public: // ======================================================================== @@ -1604,25 +1332,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** get variables from dataset in form of the table - * @param data input dataset - * @param vars list of variables - * @param cuts selection criteria - * @param table output table - * @param weights column of weigths (empty for non-weighted data) - * @param first first entry - * @param last last entry - */ - static unsigned long - get_table - ( const RooAbsData* data , - const Names& vars , - const TCut& cuts , - Table& table , - Column& weights , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== /** get variables from dataset in form of the table * @param data input dataset * @param vars list of variables @@ -1663,27 +1372,6 @@ namespace Ostap const unsigned long first = 0 , const unsigned long last = LAST ) ; // ======================================================================== - /** get variables from dataset in form of the table - * @param data input dataset - * @param vars list of variables - * @param cuts selection criteria - * @param table output table - * @param weights column of weigths (empty for non-weighted data) - * @param cutrange cut range - * @param first first entry - * @param last last entry - */ - static unsigned long - get_table - ( const RooAbsData* data , - const Names& vars , - const TCut& cuts , - Table& table , - Column& weights , - const std::string& cutrange , - const unsigned long first = 0 , - const unsigned long last = LAST ) ; - // ======================================================================== } ; // ========================================================================== } // end of namespace Ostap diff --git a/source/include/Ostap/WStatEntity.h b/source/include/Ostap/WStatEntity.h index dd65cabb..47c81bf6 100644 --- a/source/include/Ostap/WStatEntity.h +++ b/source/include/Ostap/WStatEntity.h @@ -29,16 +29,18 @@ namespace Ostap public: // the basic getters // ====================================================================== /// total number of entries - unsigned long long n () const { return m_weights.n() ; } + unsigned long long n () const { return m_weights.n() ; } /// the first weighted moment/mean-value - double mu () const { return m_mu ; } + double mu () const { return m_mu ; } /// the second central weighted moment/dispersion/variance - double mu2 () const { return m_mu2 ; } + double mu2 () const { return m_mu2 ; } // ====================================================================== public: // derived getters and aliases // ====================================================================== + /// empty ? + bool empty () const { return m_weights.empty() ; } /// get the actual number of entries - unsigned long long nEntries () const { return n() ; } + unsigned long long nEntries () const { return n () ; } /// mean-value double mean () const { return m_mu ; } /// error im mean-value @@ -100,6 +102,8 @@ namespace Ostap // ======================================================================== /// reset the counters void reset () ; + /// swap two cunters + void swap ( WStatEntity& right ) ; /// representation as string std::string toString () const; /// printout to std::ostream @@ -171,6 +175,9 @@ namespace Ostap /// conversion to string inline std::string to_string ( const WStatEntity& e ) { return e.toString() ;} // ========================================================================== + /// swap two counters + inline void swap ( WStatEntity& a , WStatEntity& b ) { a.swap( b ) ; } + // ========================================================================= } // The end of namespace Ostap // ============================================================================ // The END diff --git a/source/src/Covariance.cpp b/source/src/Covariance.cpp new file mode 100644 index 00000000..c9eb7db3 --- /dev/null +++ b/source/src/Covariance.cpp @@ -0,0 +1,278 @@ +// ============================================================================ +// Include files +// ============================================================================ +// STD & STL +// ============================================================================ +#include +#include +// ============================================================================ +// Ostap +// ============================================================================ +#include "Ostap/Math.h" +#include "Ostap/StatEntity.h" +#include "Ostap/WStatEntity.h" +#include "Ostap/Covariance.h" +// ============================================================================ +// local +// ============================================================================ +#include "Exception.h" +// ============================================================================ +/** @file + * Implementation file for class Ostap::Math::Covariance + * @see Ostap::Math::Covariance + * @see Ostap::Math::WCovariance + * @author Vanya Belyaev Ivan.Belyaev@cern.ch + * @date 2024-07-22 + */ +// ============================================================================ +namespace +{ + // ========================================================================== + const Ostap::Math::Zero s_zero {} ; + const Ostap::Math::Equal_To s_equal {} ; // equality criteria for doubles + // ========================================================================== +} +// ============================================================================ +// constructor from two counters and second moment +// ============================================================================ +Ostap::Math::Covariance::Covariance +( const Ostap::Math::Covariance::Counter& cnt1 , + const Ostap::Math::Covariance::Counter& cnt2 , + const double corr ) + : m_cnt1 ( cnt1 ) + , m_cnt2 ( cnt2 ) + , m_cov2m ( corr ) +{ + Ostap::Assert ( cnt1.n () == cnt2.n() , + "Ostap::Math::Covariance: invalid counters!" , + "Ostap::Math::Covariance" ) ; + // + const double acorr = std::abs ( corr ) ; + Ostap::Assert ( acorr <= 1.0 || s_equal ( acorr , 1.0 ) , + "Ostap::Math::Covariance: invalid correlation!" , + "Ostap::Math::Covariance" ) ; + // + if ( 1 < corr ) { m_cov2m = 1 ; } + else if ( -1 > corr ) { m_cov2m = -1 ; } + // + if ( empty() ) { m_cov2m = 0 ; } + else + { + const double covx = m_cnt1.mu2() ; + const double covy = m_cnt2.mu2() ; + if ( s_zero ( covx ) || s_zero ( covy ) ) { m_cov2m = 0 ; } + Ostap::Assert ( 0<= covx & 0 <= covy , + "Ostap::Math::Covariance: invalid variances!" , + "Ostap::Math::Covariance" ) ; + m_cov2m *= std::sqrt ( covx * covy ) ; + m_cov2m *= n () ; + } +} +// ============================================================================ +// add two counters +// ============================================================================ +Ostap::Math::Covariance& +Ostap::Math::Covariance::add +( const Ostap::Math::Covariance& right ) +{ + if ( right.empty () ) { return *this ; } + else if ( empty () ) { *this = right ; return *this ; } + // + const unsigned long long wA = n () ; + const unsigned long long wB = right.n () ; + const unsigned long long wW = wA + wB ; + // + const double xA = m_cnt1.mean () ; + const double xB = right.m_cnt1.mean () ; + const double yA = m_cnt2.mean () ; + const double yB = right.m_cnt2.mean () ; + // + m_cov2m += right.m_cov2m + ( xB - xA ) * ( yB - yA ) * wA * wB / wW ; + // + m_cnt1 += right.m_cnt1 ; + m_cnt2 += right.m_cnt2 ; + // + return *this ; +} +// ========================================================== +double Ostap::Math::Covariance::correlation () const +{ + if ( empty() ) { return 0; } + // + const double xv = m_cnt1.mu2() ; + const double yv = m_cnt2.mu2() ; + const double cc = covariance() ; + // + if ( s_zero ( cc ) ) { return 0 ; } + if ( xv < 0 || yv < 0 ) { return 0 ; } + // + return cc / std::sqrt ( xv * yv ) ; +} +// ============================================================================ +// get the covariance matrix +// ============================================================================ +Ostap::Math::Covariance::Matrix +Ostap::Math::covariance +( const Ostap::Math::Covariance& c ) +{ + // + if ( c.empty() ) { return Ostap::Math::Covariance::Matrix() ; } + // + static std::array buffer; + buffer [ 0 ] = c.counter1().mu2() ; + buffer [ 1 ] = c.covariance () ; + buffer [ 2 ] = c.counter2().mu2() ; + return Ostap::Math::Covariance::Matrix ( buffer.begin() , buffer.end() ) ; +} +// ============================================================================ +// reset counters +// ============================================================================ +void Ostap::Math::Covariance::reset () +{ + m_cnt1.reset() ; + m_cnt2.reset() ; + m_cov2m = 0 ; +} +// ============================================================================ +// get the correlation matrix +// ============================================================================ +Ostap::Math::Covariance::Matrix +Ostap::Math::correlation +( const Ostap::Math::Covariance& c ) +{ + // + if ( c.empty() ) { return Ostap::Math::Covariance::Matrix() ; } + // + static std::array buffer; + buffer [ 0 ] = 1 ; + buffer [ 1 ] = c.correlation () ; + buffer [ 2 ] = 1 ; + return Ostap::Math::Covariance::Matrix ( buffer.begin() , buffer.end() ) ; +} +// ============================================================================ +// constructor from two counters and second moment +// ============================================================================ +Ostap::Math::WCovariance::WCovariance +( const Ostap::Math::WCovariance::Counter& cnt1 , + const Ostap::Math::WCovariance::Counter& cnt2 , + const double corr ) + : m_cnt1 ( cnt1 ) + , m_cnt2 ( cnt2 ) + , m_cov2m ( corr ) +{ + Ostap::Assert ( cnt1.weights () == cnt2.weights() , + "Ostap::Math::WCovariance: invalid counters!" , + "Ostap::Math::WCovariance" ) ; + // + const double acorr = std::abs ( corr ) ; + Ostap::Assert ( acorr <= 1.0 || s_equal ( acorr , 1.0 ) , + "Ostap::Math::Covariance: invalid correlation!" , + "Ostap::Math::Covariance" ) ; + // + if ( 1 < corr ) { m_cov2m = 1 ; } + else if ( -1 > corr ) { m_cov2m = -1 ; } + // + if ( empty() ) { m_cov2m = 0 ; } + else + { + const double covx = m_cnt1.mu2() ; + const double covy = m_cnt2.mu2() ; + if ( s_zero ( covx ) || s_zero ( covy ) ) { m_cov2m = 0 ; } + Ostap::Assert ( 0<= covx & 0 <= covy , + "Ostap::Math::Covariance: invalid variances!" , + "Ostap::Math::Covariance" ) ; + m_cov2m *= std::sqrt ( covx * covy ) ; + m_cov2m *= w () ; + } +} +// ============================================================================ +// add two counters +// ============================================================================ +Ostap::Math::WCovariance& +Ostap::Math::WCovariance::add +( const Ostap::Math::WCovariance& right ) +{ + if ( right.empty () ) { return *this ; } + else if ( empty () ) { *this = right ; return *this ; } + // + const long double wA = w () ; + const long double wB = right.w () ; + const long double wW = wA + wB ; + // + m_cnt1 += right.m_cnt1 ; + m_cnt2 += right.m_cnt2 ; + // + const double xA = m_cnt1.mean () ; + const double xB = right.m_cnt1.mean () ; + const double yA = m_cnt2.mean () ; + const double yB = right.m_cnt2.mean () ; + // + m_cov2m += right.m_cov2m + ( xB - xA ) * ( yB - yA ) * wA * wB / wW ; + // + m_cnt1 += right.m_cnt1 ; + m_cnt2 += right.m_cnt2 ; + // + return *this ; +} +// ============================================================================ +// reset counters +// ============================================================================ +void Ostap::Math::WCovariance::reset () +{ + m_cnt1.reset() ; + m_cnt2.reset() ; + m_cov2m = 0 ; +} +// ========================================================== +double Ostap::Math::WCovariance::correlation () const +{ + if ( empty() ) { return 0; } + // + const double xv = m_cnt1.mu2() ; + const double yv = m_cnt2.mu2() ; + const double cc = covariance() ; + // + if ( s_zero ( cc ) ) { return 0 ; } + if ( xv < 0 || yv < 0 ) { return 0 ; } + // + return cc / std::sqrt ( xv * yv ) ; +} +// ============================================================================ +// get the covariance matrix +// ============================================================================ +Ostap::Math::WCovariance::Matrix +Ostap::Math::covariance +( const Ostap::Math::WCovariance& c ) +{ + // + if ( c.empty() ) { return Ostap::Math::WCovariance::Matrix() ; } + // + static std::array buffer; + buffer [ 0 ] = c.counter1().mu2() ; + buffer [ 1 ] = c.covariance () ; + buffer [ 2 ] = c.counter2().mu2() ; + return Ostap::Math::WCovariance::Matrix ( buffer.begin() , buffer.end() ) ; +} +// ============================================================================ +// get the correlation matrix +// ============================================================================ +Ostap::Math::WCovariance::Matrix +Ostap::Math::correlation +( const Ostap::Math::WCovariance& c ) +{ + // + if ( c.empty() ) { return Ostap::Math::WCovariance::Matrix() ; } + // + static std::array buffer; + buffer [ 0 ] = 1 ; + buffer [ 1 ] = c.correlation () ; + buffer [ 2 ] = 1 ; + return Ostap::Math::WCovariance::Matrix ( buffer.begin() , buffer.end() ) ; +} + + + + +// ============================================================================ +// The END +// ============================================================================ diff --git a/source/src/StatEntity.cpp b/source/src/StatEntity.cpp index 6401d055..89b75066 100644 --- a/source/src/StatEntity.cpp +++ b/source/src/StatEntity.cpp @@ -219,6 +219,17 @@ void Ostap::StatEntity::reset() m_max = - std::numeric_limits::max() ; } // ============================================================================ +// swap two counters +// ============================================================================ +void Ostap::StatEntity::swap ( Ostap::StatEntity& right ) +{ + std::swap ( m_n , right.m_n ) ; + std::swap ( m_mu , right.m_mu ) ; + std::swap ( m_mu2 , right.m_mu2 ) ; + std::swap ( m_min , right.m_min ) ; + std::swap ( m_max , right.m_max ) ; +} +// ============================================================================ // representation as string // ============================================================================ std::string Ostap::StatEntity::toString () const diff --git a/source/src/StatVar.cpp b/source/src/StatVar.cpp index aadd6f9d..ba7eca2b 100644 --- a/source/src/StatVar.cpp +++ b/source/src/StatVar.cpp @@ -14,7 +14,6 @@ // ============================================================================ #include "RVersion.h" #include "TTree.h" -#include "TCut.h" #include "TMatrixTSym.h" #include "RooDataSet.h" #include "RooAbsReal.h" @@ -970,24 +969,6 @@ bool Ostap::StatVar::hasEntry return false ; } // ============================================================================ -/* check if there is at least one entry that satisfies criteria - * @param tree (INPUT) the tree - * @param cuts (INPUT) criteria - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ -// ============================================================================ -bool Ostap::StatVar::hasEntry -( TTree* tree , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return hasEntry ( tree , _cuts , first , last ) ; -} -// ============================================================================ /** check if there is at least one entry that satisfies criteria * @param data (INPUT) data * @param cuts (INPUT) criteria @@ -1044,45 +1025,6 @@ bool Ostap::StatVar::hasEntry const unsigned long last ) { return hasEntry ( data , cuts , "" , first , last ) ; } // ============================================================================ -/** check if there is at least one entry that satisfies criteria - * @param data (INPUT) data - * @param cuts (INPUT) criteria - * @param cut_range (INPUT) cut range - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ -// ============================================================================ -bool Ostap::StatVar::hasEntry -( const RooAbsData* data , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return hasEntry ( data , _cuts , cut_range , first , last ) ; -} -// ============================================================================ -/** check if there is at least one entry that satisfies criteria - * @param data (INPUT) data - * @param cuts (INPUT) criteria - * @param cut_range (INPUT) cut range - * @param first (INPUT) the first entry - * @param last (INPUT) the last entry - * @return true if there exist at leats one entry - */ -// ============================================================================ -bool Ostap::StatVar::hasEntry -( const RooAbsData* data , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return hasEntry ( data , _cuts , "" , first , last ) ; -} -// ============================================================================ /* build statistic for the expression * @param tree (INPUT) the tree * @param expression (INPUT) the expression @@ -1096,14 +1038,14 @@ bool Ostap::StatVar::hasEntry * @date 2013-10-13 */ // ============================================================================ -Ostap::StatVar::Statistic +Ostap::StatEntity Ostap::StatVar::statVar ( TTree* tree , const std::string& expression , const unsigned long first , const unsigned long last ) { - Statistic result ; + Ostap::StatEntity result {} ; if ( 0 == tree || last <= first ) { return result ; } // RETURN Ostap::Formula formula ( expression , tree ) ; if ( !formula.GetNdim() ) { return result ; } // RETURN @@ -1189,32 +1131,6 @@ Ostap::StatVar::statVar return result ; } // ============================================================================ -/* build statistic for the expression - * @param tree (INPUT) the tree - * @param expression (INPUT) the expression - * @param cuts (INPUT) the selection criteria - * - * @code - * tree = ... - * stat = tree.statVar( 'S_sw' ,'pt>1000') - * @endcode - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2013-10-13 - */ -// ============================================================================ -Ostap::StatVar::Statistic -Ostap::StatVar::statVar -( TTree* tree , - const std::string& expression , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statVar ( tree , expression , _cuts , first , last ) ; -} -// ============================================================================ /* build statistic for the expressions * @param tree (INPUT) the tree * @param result (UPDATE) the output statistics for specified expressions @@ -1356,242 +1272,155 @@ unsigned long Ostap::StatVar::statVars return results.empty() ? 0 : result[0].nEntries() ; } // ============================================================================ -/** build statistic for the expressions - * @param tree (INPUT) the tree - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - * @param cuts (INPUT) the selection criteria - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-11-04 - */ -// ============================================================================ -unsigned long Ostap::StatVar::statVars -( TTree* tree , - std::vector& result , - const Ostap::StatVar::Names& expressions , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statVars ( tree , result , expressions , _cuts , first , last ) ; -} -// ============================================================================ /* calculate the covariance of two expressions * @param tree (INPUT) the input tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 + * @param exp1 (INPUT) the first expression + * @param exp2 (INPUT) the second expression + * @return covariance + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2024-07-22 */ // ============================================================================ -unsigned long +Ostap::StatVar::Covariance Ostap::StatVar::statCov ( TTree* tree , const std::string& exp1 , const std::string& exp2 , - Ostap::StatVar::Statistic& stat1 , - Ostap::StatVar::Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , const unsigned long first , const unsigned long last ) { // - stat1.reset () ; - stat2.reset () ; - Ostap::Math::setToScalar ( cov2 , 0.0 ) ; - // - if ( 0 == tree || last <= first ) { return 0 ; } // RETURN + Ostap::Assert ( nullptr != tree , + "Invalid TTree" , + "Ostap::StatVar::statCov" ) ; + // + // prepare the result + Ostap::StatVar::Covariance result {} ; + /// + if ( last <= first ) { return result ; } // RETURN + Ostap::Formula formula1 ( exp1 , tree ) ; - if ( !formula1 .ok () ) { return 0 ; } // RETURN + Ostap::Assert ( formula1.ok() , + std::string ( "Invalid first expression: " ) + exp1 , + "Ostap::StatVar::statCov" ) ; Ostap::Formula formula2 ( exp2 , tree ) ; - if ( !formula2 .ok () ) { return 0 ; } // RETURN + Ostap::Assert ( formula2.ok() , + std::string ( "Invalid second expression: " ) + exp2 , + "Ostap::StatVar::statCov" ) ; // Ostap::Utils::Notifier notify ( tree , &formula1 , &formula2 ) ; // - const unsigned long nEntries = - std::min ( last , (unsigned long) tree->GetEntries() ) ; + const unsigned long nEntries = std::min ( last , (unsigned long) tree->GetEntries() ) ; // std::vector results1 {} ; std::vector results2 {} ; + // for ( unsigned long entry = first ; entry < nEntries ; ++entry ) - { - // - long ievent = tree->GetEntryNumber ( entry ) ; - if ( 0 > ievent ) { break ; } // BREAK - // - ievent = tree->LoadTree ( ievent ) ; - if ( 0 > ievent ) { break ; } // BREAK - // - formula1.evaluate ( results1 ) ; - formula2.evaluate ( results2 ) ; - // - for ( const long double v1 : results1 ) - { - for ( const long double v2 : results2 ) - { - // - stat1 += v1 ; - stat2 += v2 ; - // - cov2 ( 0 , 0 ) += v1*v1 ; - cov2 ( 0 , 1 ) += v1*v2 ; - cov2 ( 1 , 1 ) += v2*v2 ; - } + { + // + long ievent = tree->GetEntryNumber ( entry ) ; + if ( 0 > ievent ) { break ; } // BREAK + // + ievent = tree->LoadTree ( ievent ) ; + if ( 0 > ievent ) { break ; } // BREAK + // + formula1.evaluate ( results1 ) ; + formula2.evaluate ( results2 ) ; + // + for ( const long double v1 : results1 ) + { + for ( const long double v2 : results2 ) + { + result.add ( v1 , v2 ) ; + } + } + // } - // - } - // - if ( 0 == stat1.nEntries() ) { return 0 ; } // RETURN - // - cov2 /= stat1.nEntries () ; - // - const double v1_mean = stat1.mean() ; - const double v2_mean = stat2.mean() ; - // - cov2 ( 0 , 0 ) -= v1_mean * v1_mean ; - cov2 ( 0 , 1 ) -= v1_mean * v2_mean ; - cov2 ( 1 , 1 ) -= v2_mean * v2_mean ; // - return stat1.nEntries () ; + return result ; } // ============================================================================ -/* calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection criteria - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 +/* calculate the covariance of two expressions + * @param tree (INPUT) the inpout tree + * @param exp1 (INPUT) the first expression + * @param exp2 (INPUT) the second expression + * @param cuts (INPUT) cuts/weigths expression + * @return Covariance + * @author Vanya BELYAEV Ivan.Belyaev@cern.ch + * @date 2024-07-22 */ // ============================================================================ -unsigned long +Ostap::StatVar::WCovariance Ostap::StatVar::statCov -( TTree* tree , - const std::string& exp1 , - const std::string& exp2 , - const std::string& cuts , - Ostap::StatVar::Statistic& stat1 , - Ostap::StatVar::Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const unsigned long first , - const unsigned long last ) +( TTree* tree , + const std::string& exp1 , + const std::string& exp2 , + const std::string& cuts , + const unsigned long first , + const unsigned long last ) { - if ( cuts.empty() ) - { return statCov ( tree , exp1 , exp2 , stat1 , stat2 , cov2 , first , last ) ; } // - stat1.reset () ; - stat2.reset () ; - Ostap::Math::setToScalar ( cov2 , 0.0 ) ; + Ostap::Assert ( nullptr != tree , + "Invalid TTree" , + "Ostap::StatVar::statCov" ) ; + // + // prepare the result + Ostap::StatVar::WCovariance result {} ; + /// + if ( last <= first ) { return result ; } // - if ( 0 == tree || last <= first ) { return 0 ; } // RETURN Ostap::Formula formula1 ( exp1 , tree ) ; - if ( !formula1 .ok () ) { return 0 ; } // RETURN + Ostap::Assert ( formula1.ok() , + std::string ( "Invalid first expression: " ) + exp1 , + "Ostap::StatVar::statCov" ) ; + // Ostap::Formula formula2 ( exp2 , tree ) ; - if ( !formula2 .ok () ) { return 0 ; } // RETURN - Ostap::Formula selection ( cuts , tree ) ; - if ( !selection.ok () ) { return 0 ; } // RETURN + Ostap::Assert ( formula2.ok() , + std::string ( "Invalid second expression: " ) + exp2 , + "Ostap::StatVar::statCov" ) ; + // + const std::unique_ptr selection { !cuts.empty() ? new Ostap::Formula ( cuts , tree ) : nullptr } ; + Ostap::Assert ( !selection || selection->ok () , + std::string ( "Invalid selection/weight: " ) + cuts , + "Ostap::StatVar::statCov" ) ; + /// + Ostap::Utils::Notifier notify ( tree , &formula1 , selection.get() ) ; // - Ostap::Utils::Notifier notify ( tree , &formula1 , &formula2 ) ; + const bool with_cuts = selection && selection->ok () ; // - const unsigned long nEntries = - std::min ( last , (unsigned long) tree->GetEntries() ) ; + const unsigned long nEntries = std::min ( last , (unsigned long) tree->GetEntries() ) ; // std::vector results1 {} ; std::vector results2 {} ; for ( unsigned long entry = first ; entry < nEntries ; ++entry ) - { - // - long ievent = tree->GetEntryNumber ( entry ) ; - if ( 0 > ievent ) { break ; } // RETURN - // - ievent = tree->LoadTree ( ievent ) ; - if ( 0 > ievent ) { break ; } // RETURN - // - const double w = selection.evaluate() ; - // - if ( !w ) { continue ; } // ATTENTION - // - formula1.evaluate ( results1 ) ; - formula2.evaluate ( results2 ) ; - // - for ( const long double v1 : results1 ) - { - for ( const long double v2 : results2 ) - { - // - stat1.add ( v1 , w ) ; - stat2.add ( v2 , w ) ; - // - cov2 ( 0 , 0 ) += w*v1*v1 ; - cov2 ( 0 , 1 ) += w*v1*v2 ; - cov2 ( 1 , 1 ) += w*v2*v2 ; - } + { + // + long ievent = tree->GetEntryNumber ( entry ) ; + if ( 0 > ievent ) { break ; } // RETURN + // + ievent = tree->LoadTree ( ievent ) ; + if ( 0 > ievent ) { break ; } // RETURN + // + const double w = with_cuts ? selection->evaluate() : 1.0 ; + // + if ( !w ) { continue ; } // ATTENTION + // + formula1.evaluate ( results1 ) ; + formula2.evaluate ( results2 ) ; + // + for ( const long double v1 : results1 ) + { + for ( const long double v2 : results2 ) + { + result.add ( v1 , v2 , w ) ; + } + } } - } // - if ( 0 == stat1.nEntries() || 0 == stat1.nEff () ) { return 0 ; } - // - cov2 /= stat1.weights().sum() ; - // - const double v1_mean = stat1.mean() ; - const double v2_mean = stat2.mean() ; - // - cov2 ( 0 , 0 ) -= v1_mean * v1_mean ; - cov2 ( 0 , 1 ) -= v1_mean * v2_mean ; - cov2 ( 1 , 1 ) -= v2_mean * v2_mean ; - // - return stat1.nEntries() ; -} -// ============================================================================ -/* calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection criteria - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statCov -( TTree* tree , - const std::string& exp1 , - const std::string& exp2 , - const TCut& cuts , - Ostap::StatVar::Statistic& stat1 , - Ostap::StatVar::Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - // - return statCov ( tree , - exp1 , exp2 , _cuts , - stat1 , stat2 , cov2 , - first , last ) ; + return result ; } // ============================================================================ -/* calculate the covariance of two expressions +/* calculate the covariance of several expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions * @param cuts (INPUT) the selection criteria @@ -1709,31 +1538,6 @@ Ostap::StatVar::statCov return stats[0].nEntries() ; } // ============================================================================ -/* calculate the covariance of two expressions - * @param tree (INPUT) the inpout tree - * @param vars (INPUT) expressions - * @param cuts (INPUT) the selection criteria - * @param stats (UPDATE) the statistics - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2023-02-28 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statCov -( TTree* tree , - const std::vector& vars , - const TCut& cuts , - std::vector& stats , - TMatrixTSym& cov2 , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statCov ( tree , vars , _cuts , stats , cov2 , first , last ) ; -} -// ============================================================================ /* calculate the covariance of two expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions @@ -1759,31 +1563,6 @@ Ostap::StatVar::statCov // ============================================================================ Ostap::StatVar::Statistic Ostap::StatVar::statVar -( const RooAbsData* data , - const std::string& expression , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statVar ( data , expression , _cuts , cut_range , first , last ) ; -} -// ============================================================================ -Ostap::StatVar::Statistic -Ostap::StatVar::statVar -( const RooAbsData* data , - const std::string& expression , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statVar ( data , expression , _cuts , "" , first , last ) ; -} -// ============================================================================ -Ostap::StatVar::Statistic -Ostap::StatVar::statVar ( const RooAbsData* data , const std::string& expression , const std::string& cuts , @@ -1828,58 +1607,6 @@ Ostap::StatVar::statVar return result ; } // ============================================================================ -/* build statistic for the expressions - * @param data (INPUT) input data - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - * @param cuts (INPUT) the selection - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2021-06-04 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statVars -( const RooAbsData* data , - std::vector& result , - const Ostap::StatVar::Names& expressions , - const TCut& cuts , - const unsigned long first , - const unsigned long last ) -{ - const std::string cuts_ = cuts.GetTitle() ; - return statVars ( data , result , expressions , cuts_ , "" , first , last ) ; -} -// ============================================================================ -/* build statistic for the expressions - * @param data (INPUT) input data - * @param result (UPDATE) the output statistics for specified expressions - * @param expressions (INPUT) the list of expressions - * @param cuts (INPUT) the selection - * @param cut_range (INPUT) cut range - * @param first (INPUT) the first entry to process - * @param last (INPUT) the last entry to process (not including!) - * @return number of processed entries - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2021-06-04 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statVars -( const RooAbsData* data , - std::vector& result , - const Ostap::StatVar::Names& expressions , - const TCut& cuts , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - const std::string cuts_ = cuts.GetTitle() ; - return statVars ( data , result , expressions , cuts_ , cut_range , first , last ) ; -} -// ============================================================================ /** build statistic for the expressions * @param data (INPUT) input data * @param result (UPDATE) the output statistics for specified expressions @@ -1959,92 +1686,6 @@ Ostap::StatVar::statVars return result.empty() ? 0 : result[0].nEntries() ; } // ============================================================================ -/* calculate the covariance of two expressions - * @param tree (INPUT) the input tree - * @param exp1 (INPUT) the first expresiion - * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) selection - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2014-03-27 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statCov -( const RooAbsData* data , - const std::string& exp1 , - const std::string& exp2 , - const std::string& cuts , - Ostap::StatVar::Statistic& stat1 , - Ostap::StatVar::Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - // - stat1.reset () ; - stat2.reset () ; - Ostap::Math::setToScalar ( cov2 , 0.0 ) ; - // - if ( 0 == data || last <= first ) { return 0 ; } // RETURN - // - const std::unique_ptr formula1 { make_formula ( exp1 , *data ) } ; - const std::unique_ptr formula2 { make_formula ( exp2 , *data ) } ; - const std::unique_ptr selection { make_formula ( cuts , *data , true ) } ; - // - const bool weighted = data->isWeighted() ; - const char* cutrange = cut_range.empty() ? nullptr : cut_range.c_str() ; - // - const unsigned long nEntries = std::min ( last , (unsigned long) data->numEntries() ) ; - // - for ( unsigned long entry = first ; entry < nEntries ; ++entry ) - { - // - const RooArgSet* vars = data->get( entry ) ; - if ( nullptr == vars ) { break ; } // RETURN - if ( cutrange && !vars->allInRange ( cutrange ) ) { continue ; } // CONTINUE - // - // apply cuts: - const long double wc = selection ? selection -> getVal() : 1.0L ; - if ( !wc ) { continue ; } // CONTINUE - // apply weight: - const long double wd = weighted ? data->weight() : 1.0L ; - if ( !wd ) { continue ; } // CONTINUE - // cuts & weight: - const long double w = wd * wc ; - if ( !w ) { continue ; } // CONTINUE - // - const double v1 = formula1->getVal() ; - const double v2 = formula2->getVal() ; - // - stat1.add ( v1 , w ) ; - stat2.add ( v2 , w ) ; - // - cov2 ( 0 , 0 ) += w * v1 * v1 ; - cov2 ( 0 , 1 ) += w * v1 * v2 ; - cov2 ( 1 , 1 ) += w * v2 * v2 ; - // - } - // - if ( 0 == stat1.nEntries() || 0 == stat1.nEff () ) { return 0 ; } - // - cov2 /= stat1.weights().sum() ; - // - const double v1_mean = stat1.mean() ; - const double v2_mean = stat2.mean() ; - // - cov2 ( 0 , 0 ) -= v1_mean * v1_mean ; - cov2 ( 0 , 1 ) -= v1_mean * v2_mean ; - cov2 ( 1 , 1 ) -= v2_mean * v2_mean ; - // - return stat1.nEntries() ; -} -// ============================================================================ /* calculate the covariance of several expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions @@ -2143,33 +1784,6 @@ Ostap::StatVar::statCov return stats[0].nEntries() ; } // ============================================================================ -/* calculate the covariance of several expressions - * @param tree (INPUT) the inpout tree - * @param vars (INPUT) expressions - * @param cuts (INPUT) the selection criteria - * @param stats (UPDATE) the statistics - * @param cov2 (UPDATE) the covariance matrix - * @param cut_range (INPUT) range - * @return number of processed events - * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2023-02-28 - */ -// ============================================================================ -unsigned long -Ostap::StatVar::statCov -( const RooAbsData* data , - const std::vector& vars , - const TCut& cuts , - std::vector& stats , - TMatrixTSym& cov2 , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return statCov ( data , vars , _cuts , stats , cov2 , cut_range , first , last ) ; -} -// ============================================================================ /* calculate the covariance of several expressions * @param tree (INPUT) the inpout tree * @param vars (INPUT) expressions @@ -3581,48 +3195,61 @@ Ostap::StatVar::statVar * @param frame (INPUT) data frame * @param exp1 (INPUT) the first expresiion * @param exp2 (INPUT) the second expresiion - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * + * @return covariance * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-06-18 + * @date 2024-07-22 */ // ============================================================================ -unsigned long Ostap::StatVar::statCov -( Ostap::FrameNode frame , - const std::string& exp1 , - const std::string& exp2 , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 ) -{ return statCov ( frame , exp1 , exp2 , "" , stat1 , stat2 , cov2 ) ; } +Ostap::StatVar::Covariance +Ostap::StatVar::statCov +( Ostap::FrameNode frame , + const std::string& exp1 , + const std::string& exp2 ) +{ + /// define the temporary columns + const std::string var1 = Ostap::tmp_name ( "v_" , exp1 ) ; + const std::string var2 = Ostap::tmp_name ( "v_" , exp2 ) ; + auto t = frame + .Define ( var1 , "1.0*(" + exp1 + ")" ) + .Define ( var2 , "1.0*(" + exp2 + ")" ) ; + /// + const unsigned int nSlots = std::max ( 1u , +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,22,0) + ROOT::GetThreadPoolSize () +#else + ROOT::GetImplicitMTPoolSize () +#endif + ); + // + std::vector _covs ( nSlots ? nSlots : 1 ) ; + // + auto fun = [&_covs,nSlots] ( unsigned int slot , double v1 , double v2 ) + { _covs [ slot % nSlots ].add ( v1 , v2 ) ; } ; + t.ForeachSlot ( fun , { var1 , var2 } ); + // + Covariance result ; + for ( const auto& s : _covs ) { result += s ; } + // + return result ; +} // ============================================================================ /* calculate the covariance of two expressions * @param frame (INPUT) data frame * @param exp1 (INPUT) the first expresiion * @param exp2 (INPUT) the second expresiion - * @param cuts (INPUT) the selection - * @param stat1 (UPDATE) the statistic for the first expression - * @param stat2 (UPDATE) the statistic for the second expression - * @param cov2 (UPDATE) the covariance matrix - * @return number of processed events - * + * @param cuts (INPUT) the selection/weight + * @return covariance * @author Vanya BELYAEV Ivan.Belyaev@itep.ru - * @date 2018-06-18 + * @date 2024-07-22 */ // ============================================================================ -unsigned long Ostap::StatVar::statCov +Ostap::StatVar::WCovariance +Ostap::StatVar::statCov ( Ostap::FrameNode frame , const std::string& exp1 , const std::string& exp2 , - const std::string& cuts , - Statistic& stat1 , - Statistic& stat2 , - Ostap::SymMatrix2x2& cov2 ) + const std::string& cuts ) { - // const bool no_cuts = trivial ( cuts ) ; /// define the temporary columns const std::string var1 = Ostap::tmp_name ( "v_" , exp1 ) ; @@ -3644,39 +3271,16 @@ unsigned long Ostap::StatVar::statCov #endif ); // - std::vector _sta1 ( nSlots ? nSlots : 1 ) ; - std::vector _sta2 ( nSlots ? nSlots : 1 ) ; - std::vector _cov2 ( nSlots ? nSlots : 1 ) ; + std::vector _covs ( nSlots ? nSlots : 1 ) ; // - auto fun = [&_sta1,&_sta2,&_cov2,nSlots] - ( unsigned int slot , double v1 , double v2 , double w ) { - if ( w ) - { - _sta1 [ slot % nSlots ].add ( v1 , w ) ; - _sta2 [ slot % nSlots ].add ( v2 , w ) ; - _cov2 [ slot % nSlots ]( 0 , 0 ) += w * v1 * v1 ; - _cov2 [ slot % nSlots ]( 0 , 1 ) += w * v1 * v2 ; - _cov2 [ slot % nSlots ]( 1 , 1 ) += w * v2 * v2 ; - } - } ; + auto fun = [&_covs,nSlots] ( unsigned int slot , double v1 , double v2 , double w ) + { if ( w ) { _covs [ slot % nSlots ].add ( v1 , v2 , w ) ; } } ; t.ForeachSlot ( fun , { var1 , var2 , weight } ); - // - for ( const auto& s : _sta1 ) { stat1 += s ; } - for ( const auto& s : _sta2 ) { stat2 += s ; } - for ( const auto& s : _cov2 ) { cov2 += s ; } - // - if ( 0 == stat1.nEntries() ) { return 0 ; } - // - cov2 /= stat1.nEntries() ; - // - const double v1_mean = stat1.mean() ; - const double v2_mean = stat2.mean() ; // - cov2 ( 0 , 0 ) -= v1_mean * v1_mean ; - cov2 ( 0 , 1 ) -= v1_mean * v2_mean ; - cov2 ( 1 , 1 ) -= v2_mean * v2_mean ; + WCovariance result ; + for ( const auto& s : _covs ) { result += s ; } // - return stat1.nEntries () ; + return result ; } // ============================================================================ /* calculate the moment of order "order" relative to the center "center" @@ -4372,30 +3976,6 @@ Ostap::StatVar::get_table const unsigned long last ) { return get_table ( data , vars , cuts , table , weights , "" , first , last ) ; } // ============================================================================ -/* get variables from dataset in form of the table - * @param data input dataset - * @param vars list of variables - * @param cuts selection criteria - * @param table output table - * @param weights column of weigths (empty for non-weighted data) - * @param first first entry - * @param last last entry - */ -// ============================================================================ -unsigned long -Ostap::StatVar::get_table -( const RooAbsData* data , - const Ostap::StatVar::Names& vars , - const TCut& cuts , - Ostap::StatVar::Table& table , - Ostap::StatVar::Column& weights , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return get_table ( data , vars , _cuts , table , weights , "" , first , last ) ; -} -// ============================================================================ /* get variables from dataset in form of the table * @param data input dataset * @param vars list of variables @@ -4415,31 +3995,6 @@ Ostap::StatVar::get_table const unsigned long last ) { return get_table ( data , vars , std::string("") , table , weights , "" , first , last ) ; } // ============================================================================ -/** get variables from dataset in form of the table - * @param data input dataset - * @param vars list of variables - * @param cut_range cut range - * @param table output table - * @param weights column of weigths (empty for non-weighted data) - * @param first first entry - * @param last last entry - */ -// ============================================================================ -unsigned long -Ostap::StatVar::get_table -( const RooAbsData* data , - const Ostap::StatVar::Names& vars , - const TCut& cuts , - Ostap::StatVar::Table& table , - Ostap::StatVar::Column& weights , - const std::string& cut_range , - const unsigned long first , - const unsigned long last ) -{ - const std::string _cuts = cuts.GetTitle() ; - return get_table ( data , vars , _cuts , table , weights , cut_range , first , last ) ; -} -// ============================================================================= /** get variables from dataset in form of the table * @param data input dataset * @param vars list of variables diff --git a/source/src/WStatEntity.cpp b/source/src/WStatEntity.cpp index 0bc45d41..d269c411 100644 --- a/source/src/WStatEntity.cpp +++ b/source/src/WStatEntity.cpp @@ -199,5 +199,15 @@ bool Ostap::WStatEntity::operator==( const Ostap::WStatEntity& right ) const std::tie ( right.m_mu , right.m_mu2 , right.m_values, right.m_weights ) ; } // ============================================================================ +// swap two counters +// ============================================================================ +void Ostap::WStatEntity::swap ( Ostap::WStatEntity& right ) +{ + std::swap ( m_mu , right.m_mu ) ; + std::swap ( m_mu2 , right.m_mu2 ) ; + Ostap::swap ( m_values , right.m_values ) ; + Ostap::swap ( m_weights , right.m_weights ) ; +} +// ============================================================================ // The END // ============================================================================ diff --git a/source/src/dict/Ostap.hh b/source/src/dict/Ostap.hh index 466c42b5..b6b80731 100644 --- a/source/src/dict/Ostap.hh +++ b/source/src/dict/Ostap.hh @@ -23,6 +23,7 @@ #include "Ostap/Choose.h" #include "Ostap/Clenshaw.h" #include "Ostap/Combine.h" +#include "Ostap/Covariance.h" #include "Ostap/Dalitz.h" #include "Ostap/DalitzIntegrator.h" #include "Ostap/DataFrameActions.h"