00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <TH1D.h>
00011 #include <TF1.h>
00012 #include <TGraph.h>
00013 #include <TString.h>
00014 #include <TPad.h>
00015 #include <TRandom3.h>
00016 #include <TLegend.h>
00017
00018 #include "BAT/BCLog.h"
00019 #include "BAT/BCDataSet.h"
00020 #include "BAT/BCDataPoint.h"
00021 #include "BAT/BCMath.h"
00022
00023 #include "BAT/BCHistogramFitter.h"
00024
00025
00026
00027 BCHistogramFitter::BCHistogramFitter() : BCModel("HistogramFitter")
00028 {
00029 fHistogram = 0;
00030 fFitFunction = 0;
00031
00032
00033 this -> MCMCSetNIterationsRun(2000);
00034 this -> SetFillErrorBand();
00035 fFlagIntegration = true;
00036 }
00037
00038
00039
00040 BCHistogramFitter::BCHistogramFitter(TH1D * hist, TF1 * func) : BCModel("HistogramFitter")
00041 {
00042 fHistogram = 0;
00043 fFitFunction = 0;
00044
00045 this -> SetHistogram(hist);
00046 this -> SetFitFunction(func);
00047
00048 this -> MCMCSetNIterationsRun(2000);
00049 this -> SetFillErrorBand();
00050 fFlagIntegration = true;
00051 }
00052
00053
00054
00055 int BCHistogramFitter::SetHistogram(TH1D * hist)
00056 {
00057
00058 if (!hist)
00059 {
00060 BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetHistogram() : TH1D not created.");
00061 return 0;
00062 }
00063
00064
00065 fHistogram = hist;
00066
00067
00068
00069
00070 BCDataSet * ds = new BCDataSet();
00071
00072
00073 int nbins = fHistogram -> GetNbinsX();
00074 for (int i = 0; i < nbins; ++i)
00075 {
00076 BCDataPoint* dp = new BCDataPoint(2);
00077 ds -> AddDataPoint(dp);
00078 }
00079
00080
00081 this -> SetDataSet(ds);
00082
00083
00084 double xmin = hist -> GetBinLowEdge(1);
00085 double xmax = hist -> GetBinLowEdge(nbins+1);
00086
00087
00088 double histymin = hist -> GetMinimum();
00089 double histymax = hist -> GetMaximum();
00090
00091
00092
00093 double ymin = TMath::Max(0., histymin - 5.*sqrt(histymin));
00094 double ymax = histymax + 5.*sqrt(histymax);
00095
00096
00097 this -> SetDataBoundaries(0, xmin, xmax);
00098 this -> SetDataBoundaries(1, ymin, ymax);
00099
00100
00101 this -> SetFitFunctionIndices(0, 1);
00102
00103
00104 return 1;
00105 }
00106
00107
00108
00109 int BCHistogramFitter::SetFitFunction(TF1 * func)
00110 {
00111
00112 if(!func)
00113 {
00114 BCLog::Out(BCLog::error,BCLog::error,"BCHistogramFitter::SetFitFunction() : TF1 not created.");
00115 return 0;
00116 }
00117
00118
00119 fFitFunction = func;
00120
00121
00122 this -> SetName(TString::Format("HistogramFitter with %s",fFitFunction->GetName()));
00123
00124
00125 fParameterSet -> clear();
00126
00127
00128 int n = func -> GetNpar();
00129
00130
00131 for (int i = 0; i < n; ++i)
00132 {
00133 double xmin;
00134 double xmax;
00135
00136 func -> GetParLimits(i, xmin, xmax);
00137
00138 this -> AddParameter(func->GetParName(i), xmin, xmax);
00139 }
00140
00141
00142 return 1;
00143 }
00144
00145
00146
00147 BCHistogramFitter::~BCHistogramFitter()
00148 {}
00149
00150
00151
00152 double BCHistogramFitter::LogAPrioriProbability(std::vector <double> parameters)
00153 {
00154
00155 double logprob = 0.;
00156 for(unsigned int i=0; i < this -> GetNParameters(); i++)
00157 logprob -= log(this -> GetParameter(i) -> GetRangeWidth());
00158
00159 return logprob;
00160 }
00161
00162
00163
00164 double BCHistogramFitter::LogLikelihood(std::vector <double> params)
00165 {
00166
00167 double loglikelihood = 0;
00168
00169
00170 fFitFunction -> SetParameters(¶ms[0]);
00171
00172
00173 int nbins = fHistogram -> GetNbinsX();
00174
00175
00176 double dx = fHistogram -> GetBinWidth(1);
00177
00178
00179 double fedgelow = fFitFunction -> Eval(fHistogram -> GetBinLowEdge(1));
00180
00181
00182 for (int ibin = 1; ibin <= nbins; ++ibin)
00183 {
00184
00185 double xedgehi = fHistogram -> GetBinLowEdge(ibin) + dx;
00186
00187
00188 double fedgehi = fFitFunction -> Eval(xedgehi);
00189
00190
00191 double y = fHistogram -> GetBinContent(ibin);
00192
00193 double yexp = 0.;
00194
00195
00196 if (fFlagIntegration)
00197 yexp = fFitFunction -> Integral(xedgehi-dx, xedgehi);
00198
00199
00200 else
00201 {
00202 yexp = fedgelow * dx + (fedgehi - fedgelow) * dx / 2.;
00203
00204
00205 fedgelow = fedgehi;
00206 }
00207
00208
00209 loglikelihood += BCMath::LogPoisson(y, yexp);
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 return loglikelihood;
00226 }
00227
00228
00229
00230 double BCHistogramFitter::FitFunction(std::vector <double> x, std::vector <double> params)
00231 {
00232
00233 fFitFunction -> SetParameters(¶ms[0]);
00234
00235 return fFitFunction -> Eval(x[0]) * fHistogram -> GetBinWidth( fHistogram -> FindBin(x[0]) );
00236 }
00237
00238
00239
00240 int BCHistogramFitter::Fit(TH1D * hist, TF1 * func)
00241 {
00242
00243 if (hist)
00244 this -> SetHistogram(hist);
00245 else
00246 {
00247 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Histogram not defined.");
00248 return 0;
00249 }
00250
00251
00252 if (func)
00253 this -> SetFitFunction(func);
00254 else
00255 {
00256 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Fit function not defined.");
00257 return 0;
00258 }
00259
00260
00261 this -> MarginalizeAll();
00262
00263
00264
00265 this -> FindModeMinuit( this -> GetBestFitParameters() , -1);
00266
00267
00268 double pvalue;
00269 if (this -> CalculatePValueFast(this -> GetBestFitParameters(), pvalue))
00270 fPValue = pvalue;
00271 else
00272 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::Fit() : Could not use the fast p-value evaluation.");
00273
00274
00275 this -> PrintShortFitSummary();
00276
00277
00278 return 1;
00279 }
00280
00281
00282
00283 void BCHistogramFitter::DrawFit(const char * options, bool flaglegend)
00284 {
00285 if (!fHistogram)
00286 {
00287 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Histogram not defined.");
00288 return;
00289 }
00290
00291 if (!fFitFunction)
00292 {
00293 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::DrawFit() : Fit function not defined.");
00294 return;
00295 }
00296
00297
00298 TString opt = options;
00299 opt.ToLower();
00300
00301
00302 if(!opt.Contains("same"))
00303 fHistogram -> Draw(opt.Data());
00304
00305
00306 fErrorBand = this -> GetErrorBandGraph(0.16, 0.84);
00307 fErrorBand -> Draw("f same");
00308
00309
00310 fHistogram -> Draw(TString::Format("%ssame",opt.Data()));
00311
00312
00313 fGraphFitFunction = this -> GetFitFunctionGraph( this->GetBestFitParameters() );
00314 fGraphFitFunction -> SetLineColor(kRed);
00315 fGraphFitFunction -> SetLineWidth(2);
00316 fGraphFitFunction -> Draw("l same");
00317
00318
00319 if (flaglegend)
00320 {
00321 TLegend * legend = new TLegend(0.25, 0.75, 0.55, 0.95);
00322 legend -> SetBorderSize(0);
00323 legend -> SetFillColor(kWhite);
00324 legend -> AddEntry(fHistogram, "Data", "L");
00325 legend -> AddEntry(fGraphFitFunction, "Best fit", "L");
00326 legend -> AddEntry(fErrorBand, "Error band", "F");
00327 legend -> Draw();
00328 }
00329
00330 gPad -> RedrawAxis();
00331 }
00332
00333
00334 int BCHistogramFitter::CalculatePValueFast(std::vector<double> par, double &pvalue)
00335 {
00336
00337 if (par.size() != this -> GetNParameters())
00338 {
00339 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Number of parameters is inconsistent.");
00340 return 0;
00341 }
00342
00343
00344 if (!fHistogram)
00345 {
00346 BCLog::Out(BCLog::warning, BCLog::warning,"BCHistogramFitter::CalculatePValueFast() : Histogram not defined.");
00347 return 0;
00348 }
00349
00350
00351 int nbins = fHistogram -> GetNbinsX();
00352
00353 std::vector <int> histogram;
00354 std::vector <double> expectation;
00355 histogram.assign(nbins, 0);
00356 expectation.assign(nbins, 0);
00357
00358 double logp = 0;
00359 double logp_start = 0;
00360 int counter_pvalue = 0;
00361
00362
00363 for (int ibin = 0; ibin < nbins; ++ibin)
00364 {
00365
00366 double xmin = fHistogram -> GetBinLowEdge(ibin+1);
00367 double xmax = fHistogram -> GetBinLowEdge(ibin+2);
00368
00369
00370 double yexp = fFitFunction -> Integral(xmin, xmax);
00371
00372 histogram[ibin] = int(yexp);
00373 expectation[ibin] = yexp;
00374
00375
00376 logp += BCMath::LogPoisson(double(int(yexp)), yexp);
00377 logp_start += BCMath::LogPoisson(fHistogram -> GetBinContent(ibin+1), yexp);
00378 }
00379
00380 int niter = 100000;
00381
00382
00383 for (int iiter = 0; iiter < niter; ++iiter)
00384 {
00385
00386 for (int ibin = 0; ibin < nbins; ++ibin)
00387 {
00388
00389 double ptest = fRandom -> Rndm() - 0.5;
00390
00391
00392 if (ptest > 0)
00393 {
00394
00395 double r = expectation.at(ibin) / double(histogram.at(ibin) + 1);
00396
00397
00398 if (fRandom -> Rndm() < r)
00399 {
00400 histogram[ibin] = histogram.at(ibin) + 1;
00401 logp += log(r);
00402 }
00403 }
00404
00405
00406 else if (ptest <= 0 && histogram[ibin] > 0)
00407 {
00408
00409 double r = double(histogram.at(ibin)) / expectation.at(ibin);
00410
00411
00412 if (fRandom -> Rndm() < r)
00413 {
00414 histogram[ibin] = histogram.at(ibin) - 1;
00415 logp += log(r);
00416 }
00417 }
00418
00419
00420
00421 }
00422
00423
00424 if (logp < logp_start)
00425 counter_pvalue++;
00426
00427 }
00428
00429
00430 pvalue = double(counter_pvalue) / double(niter);
00431
00432
00433 return 1;
00434 }
00435
00436