| 1 |
#include <iostream>
|
| 2 |
#include <math.h>
|
| 3 |
#include <TH1.h>
|
| 4 |
#include <TF1.h>
|
| 5 |
#include <TMatrixD.h>
|
| 6 |
#include <TMinuit.h>
|
| 7 |
|
| 8 |
|
| 9 |
|
| 10 |
// Create uncertainty band (histogram) for a given function and error matrix
|
| 11 |
// in the range of the function.
|
| 12 |
TH1D* getBand(TF1 *f, TMatrixD const& m, std::string name, bool getRelativeBand, int npx) {
|
| 13 |
|
| 14 |
double xmin = f->GetXmin();
|
| 15 |
double xmax = f->GetXmax()*1.1; //fixes problem in drawing with c option
|
| 16 |
int npar = f->GetNpar();
|
| 17 |
|
| 18 |
// Create binning (linear or log)
|
| 19 |
Double_t xvec[npx];
|
| 20 |
xvec[0] = xmin;
|
| 21 |
double dx = (xmax-xmin)/npx;
|
| 22 |
for (int i = 0; i != npx; ++i) {
|
| 23 |
xvec[i+1] = xvec[i]+dx;
|
| 24 |
}
|
| 25 |
|
| 26 |
|
| 27 |
//
|
| 28 |
// Compute partial derivatives numerically
|
| 29 |
// can be used with any fit function
|
| 30 |
//
|
| 31 |
Double_t sigmaf[npx];
|
| 32 |
TH1D* h1_band = new TH1D(name.c_str(), "", npx, xvec);
|
| 33 |
|
| 34 |
for( unsigned ipx=0; ipx<=npx; ++ipx ) {
|
| 35 |
|
| 36 |
sigmaf[ipx] = 0.;
|
| 37 |
Double_t partDeriv[npar];
|
| 38 |
|
| 39 |
//compute partial derivatives of f wrt its parameters:
|
| 40 |
for( unsigned ipar=0; ipar<npar; ++ipar ) {
|
| 41 |
|
| 42 |
Float_t pi = f->GetParameter(ipar);
|
| 43 |
Float_t dpi = sqrt(m[ipar][ipar])*0.01; //small compared to the par sigma
|
| 44 |
f->SetParameter(ipar, pi+dpi);
|
| 45 |
Float_t fplus = f->Eval(xvec[ipx]-0.5*dx);
|
| 46 |
f->SetParameter(ipar, pi-dpi);
|
| 47 |
Float_t fminus = f->Eval(xvec[ipx]-0.5*dx);
|
| 48 |
f->SetParameter(ipar, pi); //put it back as it was
|
| 49 |
|
| 50 |
partDeriv[ipar] = (fplus-fminus)/(2.*dpi);
|
| 51 |
|
| 52 |
} //for params
|
| 53 |
|
| 54 |
//compute sigma(f) at x:
|
| 55 |
for( unsigned ipar=0; ipar<npar; ++ipar ) {
|
| 56 |
for( unsigned jpar=0; jpar<npar; ++jpar ) {
|
| 57 |
sigmaf[ipx] += partDeriv[ipar]*partDeriv[jpar]*m[ipar][jpar];
|
| 58 |
}
|
| 59 |
}
|
| 60 |
sigmaf[ipx] = sqrt(sigmaf[ipx]); //absolute band
|
| 61 |
|
| 62 |
h1_band->SetBinContent( ipx, f->Eval(xvec[ipx]-0.5*dx) );
|
| 63 |
if( getRelativeBand )
|
| 64 |
h1_band->SetBinError( ipx, sigmaf[ipx]/f->Eval(xvec[ipx]-0.5*dx) );
|
| 65 |
else
|
| 66 |
h1_band->SetBinError( ipx, sigmaf[ipx] );
|
| 67 |
|
| 68 |
} //for points
|
| 69 |
|
| 70 |
|
| 71 |
return h1_band;
|
| 72 |
|
| 73 |
} //getband
|
| 74 |
|
| 75 |
TH1D* getBand( TF1* f, const std::string& name ) {
|
| 76 |
const int ndim_resp_q = f->GetNpar();
|
| 77 |
TMatrixD emat_resp_q(ndim_resp_q, ndim_resp_q);
|
| 78 |
gMinuit->mnemat(&emat_resp_q[0][0], ndim_resp_q);
|
| 79 |
|
| 80 |
return getBand(f, emat_resp_q, name,false,1000);
|
| 81 |
}
|
| 82 |
|
