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() ;

}