library: libHist #include "TFractionFitter.h"

# TFractionFitter

class description - source file - inheritance tree (.pdf)

## class TFractionFitter : public TObject

| Inheritance Chart:
`TObject`
| <- | `TFractionFitter` |

        private:

          void CheckConsistency()
          void CheckParNo(Int_t parm) const
          void ComputeChisquareLambda()
          void ComputeFCN(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t flag)
          void FindPrediction(int bin, double* fractions, double& Ti, int& k0, double& Aki) const
          void GetRanges(Int_t& minX, Int_t& maxX, Int_t& minY, Int_t& maxY, Int_t& minZ, Int_t& maxZ) const

        public:

                          TFractionFitter()
                          TFractionFitter(TH1* data, TObjArray* MCs)
                  virtual ~TFractionFitter()
           static TClass* Class()
                     void Constrain(Int_t parm, Double_t low, Double_t high)
                     void ErrorAnalysis(Double_t UP)
                    Int_t Fit()
                 Double_t GetChisquare() const
          TVirtualFitter* GetFitter() const
                     TH1* GetMCPrediction(Int_t parm) const
                    Int_t GetNDF() const
                     TH1* GetPlot()
                 Double_t GetProb() const
                     void GetResult(Int_t parm, Double_t& value, Double_t& error) const
          virtual TClass* IsA() const
                     void ReleaseRangeX()
                     void ReleaseRangeY()
                     void ReleaseRangeZ()
                     void SetData(TH1* data)
                     void SetMC(Int_t parm, TH1* MC)
                     void SetRangeX(Int_t low, Int_t high)
                     void SetRangeY(Int_t low, Int_t high)
                     void SetRangeZ(Int_t low, Int_t high)
                     void SetWeight(Int_t parm, TH1* weight)
             virtual void ShowMembers(TMemberInspector& insp, char* parent)
             virtual void Streamer(TBuffer& b)
                     void StreamerNVirtual(TBuffer& b)
                     void UnConstrain(Int_t parm)

### Data Members

        protected:

             Bool_t fFitDone       flags whether a valid fit has been performed
              Int_t fLowLimitX     first bin in X dimension
              Int_t fHighLimitX    last  bin in X dimension
              Int_t fLowLimitY     first bin in Y dimension
              Int_t fHighLimitY    last  bin in Y dimension
              Int_t fLowLimitZ     first bin in Z dimension
              Int_t fHighLimitZ    last  bin in Z dimension
              Int_t fNpfits        Number of points used in the fit
              Int_t fNDF           Number of degrees of freedom in the fit
           Double_t fChisquare     Template fit chisquare
          TObjArray fAji           array of pointers to predictions of real template distributions
               TH1* fData          pointer to the "data" histogram to be fitted to
          TObjArray fMCs           array of pointers to template histograms
          TObjArray fWeights       array of pointers to corresponding weight factors (may be null)
           Double_t fIntegralData  "data" histogram content integral over allowed fit range
          Double_t* fIntegralMCs   same for template histograms (weights not taken into account)
          Double_t* fFractions     template fractions scaled to the "data" histogram statistics
               TH1* fPlot          pointer to histogram containing summed template predictions
              Int_t fNpar          number of fit parameters

## Class Description

     Fits MC fractions to data histogram (a la HMCMLL, see R. Barlow and C. Beeston,
     Comp. Phys. Comm. 77 (1993) 219-228, and http://www.hep.man.ac.uk/~roger/hfrac.f).

     The virtue of this fit is that it takes into account both data and Monte Carlo
     statistical uncertainties. The way in which this is done is through a standard
     likelihood fit using Poisson statistics; however, the template (MC) predictions
     are also varied within statistics, leading to additional contributions to the
     overall likelihood. This leads to many more fit parameters (one per bin per
     template), but the minimisation with respect to these additional parameters is
     done analytically rather than introducing them as formal fit parameters. Some
     special care needs to be taken in the case of bins with zero content. For more
     details please see the original publication cited above.

     An example application of this fit is given below. For a TH1* histogram
     ("data") fitted as the sum of three Monte Carlo sources ("mc"):

     {
       TH1F *data;                              //data histogram
       TH1F *mc0;                               // first MC histogram
       TH1F *mc1;                               // second MC histogram
       TH1F *mc2;                               // third MC histogram
       ....                                     // retrieve histograms
       TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
       mc->Add(mc0);
       mc->Add(mc1);
       mc->Add(mc2);
       TFractionFitter* fit = new TFractionFitter(data, mc); // initialise
       fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
       fit->SetRangeX(1,15);                    // use only the first 15 bins in the fit
       Int_t status = fit->Fit();               // perform the fit
       cout << "fit status: " << status << endl;
       if (status == 0) {                       // check on fit status
         TH1F* prediction = (TH1F*) fit->GetPlot();
         data->Draw("Ep");
         result->Draw("same");
       }
     }

     Instantiation
     =============
     A fit object is instantiated through
         TFractionFitter* fit = new TFractionFitter(data, mc);
     A number of basic checks (intended to ensure that the template histograms
     represent the same "kind" of distribution as the data one) are carried out.
     The TVirtualFitter object is then addressed and all fit parameters (the
     template fractions) declared (initially unbounded).

     Applying constraints
     ====================
     Fit parameters can be constrained through
         fit->Constrain(parameter #, lower bound, upper bound);
     Setting lower bound = upper bound = 0 removes the constraint (a la Minuit);
     however, a function
         fit->Unconstrain(parameter #)
     is also provided to simplify this.

     Setting parameter values
     ========================
     The function
         TVirtualFitter* vFit = fit->GetFitter();
     is provided for direct access to the TVirtualFitter object. This allows to
     set and fix parameter values, and set step sizes directly.

     Restricting the fit range
     =========================
     The fit range can be restricted through
         fit->SetRangeX(first bin #, last bin #);
     and freed using
         fit->ReleaseRangeX();
     For 2D histograms the Y range can be similarly restricted using
         fit->SetRangeY(first bin #, last bin #);
         fit->ReleaseRangeY();
     and for 3D histograms also
         fit->SetRangeZ(first bin #, last bin #);
         fit->ReleaseRangeZ();

     Weights histograms
     ==================
     Weights histograms (for a motivation see the above publication) can be specified
     for the individual MC sources through
         fit->SetWeight(parameter #, pointer to weights histogram);
     and unset by specifying a null pointer.

     Obtaining fit results
     =====================
     The fit is carried out through
         Int_t status = fit->Fit();
     where  status  is the code returned from the "MINIMIZE" command. For fits
     that converged, parameter values and errors can be obtained through
         fit->GetResult(parameter #, value, error);
     and the histogram corresponding to the total Monte Carlo prediction (which
     is not the same as a simple weighted sum of the input Monte Carlo distributions)
     can be obtained by
         TH1* result = fit->GetPlot();

     Using different histograms
     ==========================
     It is possible to change the histogram being fitted through
         fit->SetData(TH1* data);
     and to change the template histogram for a given parameter number through
         fit->SetMC(parameter #, TH1* MC);
     This can speed up code in case of multiple data or template histograms;
     however, it should be done with care as any settings are taken over from
     the previous fit. In addition, neither the dimensionality nor the numbers of
     bins of the histograms should change (in that case it is better to instantiate
     a new TFractionFitter object).

     Errors
     ======
     Any serious inconsistency results in an error.

TFractionFitter() : fFitDone(kFALSE),fData(0), fPlot(0)

     TFractionFitter default constructor.

TFractionFitter(TH1* data, TObjArray *MCs) : fFitDone(kFALSE), fChisquare(0), fPlot(0)

     TFractionFitter constructor. Does a complete initialisation (including
     consistency checks, default fit range as the whole histogram but without
     under- and overflows, and declaration of the fit parameters). Note that
     the histograms are not copied, only references are used.
     Arguments:
         data: histogram to be fitted
         MCs:  array of TH1* corresponding template distributions

~TFractionFitter()

     TFractionFitter default destructor

void SetData(TH1* data)

     Change the histogram to be fitted to. Notes:
     - Parameter constraints and settings are retained from a possible previous fit.
     - Modifying the dimension or number of bins results in an error (in this case
       rather instantiate a new TFractionFitter object)

void SetMC(Int_t parm, TH1* MC)

     Change the histogram for template number <parm>. Notes:
     - Parameter constraints and settings are retained from a possible previous fit.
     - Modifying the dimension or number of bins results in an error (in this case
       rather instantiate a new TFractionFitter object)

void SetWeight(Int_t parm, TH1* weight)

     Set bin by bin weights for template number <parm> (the parameter numbering
     follows that of the input template vector).
     Weights can be "unset" by passing a null pointer.
     Consistency of the weights histogram with the data histogram is checked at
     this point, and an error in case of problems.

TVirtualFitter* GetFitter() const

     Give direct access to the underlying minimisation class. This can be
     used e.g. to modify parameter values or step sizes.

void CheckParNo(Int_t parm) const

     Function for internal use, checking parameter validity
     An invalid parameter results in an error.

void SetRangeX(Int_t low, Int_t high)

     Set the X range of the histogram to be used in the fit.
     Use ReleaseRangeX() to go back to fitting the full histogram.
     The consistency check ensures that no empty fit range occurs (and also
     recomputes the bin content integrals).
     Arguments:
         low:  lower X bin number
         high: upper X bin number

void ReleaseRangeX()

     Release restrictions on the X range of the histogram to be used in the fit.

void SetRangeY(Int_t low, Int_t high)

     Set the Y range of the histogram to be used in the fit (2D or 3D histograms only).
     Use ReleaseRangeY() to go back to fitting the full histogram.
     The consistency check ensures that no empty fit range occurs (and also
     recomputes the bin content integrals).
     Arguments:
         low:  lower Y bin number
         high: upper Y bin number

void ReleaseRangeY()

     Release restrictions on the Y range of the histogram to be used in the fit.

void SetRangeZ(Int_t low, Int_t high)

     Set the Z range of the histogram to be used in the fit (3D histograms only).
     Use ReleaseRangeY() to go back to fitting the full histogram.
     The consistency check ensures that no empty fit range occurs (and also
     recomputes the bin content integrals).
     Arguments:
         low:  lower Y bin number
         high: upper Y bin number

void ReleaseRangeZ()

     Release restrictions on the Z range of the histogram to be used in the fit.

void Constrain(Int_t parm, Double_t low, Double_t high)

     Constrain the values of parameter number <parm> (the parameter numbering
     follows that of the input template vector).
     Use UnConstrain() to remove this constraint.

void UnConstrain(Int_t parm)

     Remove the constraints on the possible values of parameter <parm>.

void CheckConsistency()

     Function used internally to check the consistency between the
     various histograms. Checks are performed on nonexistent or empty
     histograms, the precise histogram class, and the number of bins.
     In addition, integrals over the "allowed" bin ranges are computed.
     Any inconsistency results in a error.

Int_t Fit()

     Perform the fit with the default UP value.
     The value returned is the minimisation status.

void ErrorAnalysis(Double_t UP)

     Set UP to the given value (see class TMinuit), and perform a MINOS minimisation.

void GetResult(Int_t parm, Double_t& value, Double_t& error) const

     Obtain the fit result for parameter <parm> (the parameter numbering
     follows that of the input template vector).

TH1* GetPlot()

     Return the "template prediction" corresponding to the fit result (this is not
     the same as the weighted sum of template distributions, as template statistical
     uncertainties are taken into account).
     Note that the name of this histogram will simply be the same as that of the
     "data" histogram, prefixed with the string "Fraction fit to hist: ".

void GetRanges(Int_t& minX, Int_t& maxX, Int_t& minY, Int_t& maxY, Int_t& minZ, Int_t& maxZ) const

     Used internally to obtain the bin ranges according to the dimensionality of
     the histogram and the limits set by hand.

void ComputeFCN(Int_t& /*npar*/, Double_t* /*gin*/, Double_t& f, Double_t* xx, Int_t flag)

     Used internally to compute the likelihood value.

void FindPrediction(int bin, Double_t *fractions, Double_t &Ti, int& k0, Double_t &Aki) const

     Function used internally to obtain the template prediction in the individual bins

Double_t GetChisquare() const

     Return the likelihood ratio Chi-squared (chi2) for the fit.
     The value is computed when the fit is executed successfully.
     Chi2 calculation is based on the "likelihood ratio" lambda,
     lambda = L(y;n) / L(m;n),
     where L(y;n) is the likelihood of the fit result <y> describing the data <n>
     and L(m;n) is the likelihood of an unknown "true" underlying distribution
     <m> describing the data <n>. Since <m> is unknown, the data distribution is
     used instead,
     lambda = L(y;n) / L(n;n).
     Note that this ratio is 1 if the fit is perfect. The chi2 value is then
     computed according to
     chi2 = -2*ln(lambda).
     This parameter can be shown to follow a Chi-square distribution. See for
     example S. Baker and R. Cousins, "Clarification of the use of chi-square
     and likelihood functions in fits to histograms", Nucl. Instr. Meth. A221,
     pp. 437-442 (1984)

Int_t GetNDF() const

     return the number of degrees of freedom in the fit
     the fNDF parameter has been previously computed during a fit.
     The number of degrees of freedom corresponds to the number of points
     used in the fit minus the number of templates.

Double_t GetProb() const

     return the fit probability

void ComputeChisquareLambda()

     Method used internally to compute the likelihood ratio chi2
     See the function GetChisquare() for details

TH1* GetMCPrediction(Int_t parm) const

     Return the adjusted MC template (Aji) for template (parm).
     Note that the (Aji) times fractions only sum to the total prediction
     of the fit if all weights are 1.

### Inline Functions

                TClass* Class()
                TClass* IsA() const
                   void ShowMembers(TMemberInspector& insp, char* parent)
                   void Streamer(TBuffer& b)
                   void StreamerNVirtual(TBuffer& b)

Author: Frank Filthaut filthaut@hef.kun.nl 20/05/2002
Last update: root/hist:$Name: $:$Id: TFractionFitter.cxx,v 1.10 2004/08/18 12:29:44 brun Exp $