#ifndef ASIMOVHYPOTHESIS_INCLUDED_HH #define ASIMOVHYPOTHESIS_INCLUDED_HH #include #include #include #include #include #include #include #include "Math/Minimizer.h" #include "Math/Factory.h" #include "Math/Functor.h" #include "TMinuitMinimizer.h" #include "TStopwatch.h" #include "statutil.hh" #include "stringutil.hh" #include "rootutil.hh" using stringutil::print; using namespace std; /*! Class to fit simple mu* S + B models models to datasets, generate speudo-data and compute distributions of poisson-likelihood-ratio-test-statisics */ struct BinnedAnalysis2D { vector signal_histos; vector background_histos; int verb = 0; bool sparse = true; // make_pseudo_dataset( double sigstrength ) { vector r = signal_histos; // copy fill_pseudo_dataset( r, sigstrength ); return r; } /*! fill histogram with pseudo-data for the given sigstrength */ void fill_pseudo_dataset( vector& data, double sigstrength ); /*! Log of the poisson likelihood of the dataset, given the signal strength */ double log_likelihood( const vector& data , double sigstrength ); /*! The log likelihood ratio of the data as function of the tested sig. strength */ double log_likelihood_ratio( const vector& data, double sigstrength ); /*! Fit the signal-strength, optimising the log-likelihood */ double optimize_sigstrength( const vector& data , double startvalue = 1, double xmin = -3, double xmax = 1000, double stepsize = 1e-3); /*! Return the likelihood ratio test statistic for a data-set */ double test_statistic ( const vector& data ) { double s = optimize_sigstrength( data ); double r = log_likelihood_ratio( data, s ); if (verb) print( s, r ); return r; } /*! Return a histogram with the distribution of the test-statistic for the given signal strength */ TH1D test_statistic_distribution( double sig_strength , int n_pe =1000 , int nbins = 1100 , double xmin = -10 , double xmax = 100 ); /*! produce histogram to check bias of signal-strength fit */ TH2D check_signal_fit( int n_pe = 1000, double ss_min = 0, double ss_max = 10 ) { TH2D r("hsigfit","hsigfit;true signal strength;fitted-true signal strength", 50, ss_min, ss_max, 31, -5, 5 ); for (int i =0; i< n_pe; i++ ) { double ss = gRandom->Uniform( ss_min, ss_max ); auto data = make_pseudo_dataset( ss ); double s = optimize_sigstrength( data ); r.Fill( ss, (s-ss) ); } return r; } }; #endif