using namespace RooFit ;
void ExampleRangeFit()
{
// This macro demonstrates how to set up a fit in two ranges of an onservable
// in such a way that does not only fit the shapes in each region, but also
// takes into account the relative normalization of the two
//
// *** PART 1 shape fits (no extended Likelihood) ***
//
// If you perform a fit in two ranges in RooFit, e.g. pdf->fitTo(data,Range("Range1,Rang2"))
// it will construct a simple simultaneous fit of the two regions
//
// In case the pdf is not extended, i.e. a shape fit only, it will only fit the shapes in the
// two selected ranges, and not take into account the relative predicted yields in those ranges
//
// In certain models (like exponential decays) and configurations (e.g. narrow ranges that are far apart)
// the relative normalization of the ranges may carry much more information on the shape than the
// distribution inside those range, hence it is important to take that into account when doing such a fit
//
// This is particularly important for cases where the 2-range fit is meant to be representative of
// a full-range fit, but with a blinded signal region inside it.
//
//
// *** PART 2 shape+rate fits (extended likelihood) ***
//
// Also if your pdf is already extended, i.e. measuring both the distribution in the observable as well
// as the event count in the fitted region, some intervention is needed to make fits in ranges
// work in a way that corresponds to intuition
//
// If an extended fit is performed in a sub-range, the observed yield is that of the subrange, hence
// the expected event count in the model will be fitted to this smaller number and thus have a different
// interpretation. In such cases it is often preferred to keep the original interpretation of the
// observed event yield, i.e. apply a correction to the extended likelihood term in such a way
// that the interpretation of the expected event count remains that of the full range. This can
// be done by applying a correcion factor (equal to the fraction of the pdf that is contained in the
// fitted range) in the Poisson term that represents the extended likelihood term
//
// If an extended likelihood fit is performed over _two_ sub-ranges this correction factory is
// even more important: without it, each component likelihood would have a different interpretation
// of the expected event count (each corresponding to the count in its own region) and a joint
// fit of these regions with different interpretations of the same model parameter results
// in a number that is not easily interpreted.
//
// If both regions correct their interpretatin so that Nexpected is again that of the full sample
// this problem is naturally fixed and an intuitive meaning of the Nexpected is restored
//
// THE GOOD NEWS IS: if you have constructed your extended model with RooAddPdf in the
// form SumPdf = Nsig * sigPdf + Nbkg * bkgPdf THE ABOVE CORRECTION FOR THE INTEPRETATION
// OF NSIG AND NBKG IS AUTOMATICALLY APPLIED BY ROOFIT
//
RooWorkspace w("w") ;
// **** PART 1 ****
//
// Background-only fits
// Build plain exponential model
w.factory("Exponential::model(x[10,100],alpha[-0.04,-0.1,-0])") ;
// Define side band regions
w.var("x")->setRange("LEFT",10,20) ;
w.var("x")->setRange("RIGHT",60,100) ;
// Define full region
w.var("x")->setRange("FULL",10,100) ;
// Construct an extended pdf with associates a predicted event count N associated with the existing exponential model
//
// The interpretation of this event count is forcibly defined to be in the range "FULL"
// If the actual domain of x that is fitted is identical to FULL this has no affect
//
// But if the actually domain is a subset of FULL the expected event count is divided by fraction ( Int_FitRange model(x) dx / Int_FULLRange model(x) dx)
// The effect of this is that the fit will return a predicted event count in the FULL range rather than in the actually fitted range
w.factory("ExtendPdf::extmodel(model,N[0,20000],'FULL')") ;
// As a pedagogical exercise it is instructive to fit the above model to either the LEFT or RIGHT range, in both cases it will return a number compatible with
// the generated event count in the full range
//
// If we now do a simultaneous fit to the extended model, instead of the original model, the LEFT and RIGHT range will each introduce a different
// fractional correction factor for the event yield, but done is such a way that both likelihoods have a consistent definition of the meaning
// of the yield parameter N - that it is the predicted event count in the full yield
//
// This joint fit of the extmodel will now include (w.r.t the plain model fit) a product of extended terms
//
// Poisson(Nobs_LEFT | Nexp / frac_LEFT) * Poisson (Nobs_RIGHT | Nexp / frac_RIGHT )
//
// that will introduce additional sensitivity of the likelihood to the slope parameter alpha of the exponential model through the frac_LEFT and frac_RIGHT integrals
//
// In the extreme case of an exponential function and a fit in narrow LEFT and RIGHT ranges this sensitivity through the exponential term may actually be larger
// than from the shapes.
//
// This is also nicely demonstrated in the example below where the uncertainty on alpha is almost 5x smaller if the extended term is included
TCanvas* c = new TCanvas("c","c",2100,700) ;
c->Divide(3) ;
RooDataSet* d = w.pdf("model")->generate(*w.var("x"),10000) ;
c->cd(1) ;
RooFitResult* r = w.pdf("model")->fitTo(*d,Range("LEFT,RIGHT"),Save()) ;
RooPlot* frame = w.var("x")->frame() ;
d->plotOn(frame) ;
w.pdf("model")->plotOn(frame,VisualizeError(*r)) ;
w.pdf("model")->plotOn(frame) ;
w.pdf("model")->paramOn(frame,Label("B-only shape fit")) ;
frame->Draw() ;
c->cd(2) ;
RooFitResult* r2 = w.pdf("extmodel")->fitTo(*d,Range("LEFT,RIGHT"),Save()) ;
RooPlot* frame2 = w.var("x")->frame() ;
d->plotOn(frame2) ;
w.pdf("extmodel")->plotOn(frame2) ;
w.pdf("extmodel")->plotOn(frame2,VisualizeError(*r2)) ;
w.pdf("extmodel")->paramOn(frame2,Label("B-only extended fit with adjusted yield"),Layout(0.4,0.95) ) ;
frame2->Draw() ;
// *** PART 2 ***
//
// Now we repeat the above exercise, but instead of explicitly adding an extended term to a single shape pdf (RooExponential)
// we assume that we already have an extended likelihood model in the form of a RooAddPdf constructed in the form Nsig * sigPdf + Nbkg * bkgPdf
//
//
// Constructed an extended signal+background model: we add a Gaussian to the previously defined exponential background
// The signal is chosen constant in this example since it is entirely constrained in the region between the fit
// regions hence the fit in the side bands has no sensitivity to it.
c->cd(3) ;
w.factory("SUM::modelsum(Nsig[1000]*Gaussian::sig(x,mean[40],width[5]),Nbkg[10000,0,20000]*model)") ;
// This model will automatically insert the correction factor for the reinterpretation of Nsig and Nnbkg in the full ranges
//
// You should see the following lines appear on the terminal that confirm this when you fit this model
//
//[#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_LEFT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables
//[#1] INFO:Fitting -- RooAbsOptTestStatistic::ctor(nll_modelsum_modelsumData_RIGHT) fixing interpretation of coefficients of any RooAddPdf to full domain of observables
//
RooFitResult* r3 = w.pdf("modelsum")->fitTo(*d,Range("LEFT,RIGHT"),Save()) ;
RooPlot* frame3 = w.var("x")->frame() ;
d->plotOn(frame3) ;
w.pdf("modelsum")->plotOn(frame3) ;
w.pdf("modelsum")->plotOn(frame3,VisualizeError(*r3)) ;
w.pdf("modelsum")->paramOn(frame3,Label("S+B fit with auto-adjusted yields"),Layout(0.3,0.95)) ;
frame3->Draw() ;
}