Coherent sum of pdf with RooFit / RooSimultaneous

[canigia]

Hi,

is there a way to create a composite model for RooFit being a coherent sum of two functions/pdfs?

I.e. each pdf is defined as a complex-valued function, and the real valued intensity (=total pdf value) of the combined model then is the square of the sum of all complex function values for a certain input value x. The fit fraction in this case should be the integral of the intensity of each pdf alone divided by the combined intensity integral (=number of entries in the plot). The individual fit fractions then usually sum up to a different number than 1.0. This probably spoils the usual fit-fraction concept in RooFit.

However, if it is possible, I would be interested in doing even a simultaneous fit, constraining the absolute fit fraction of the signal pdf (integral of the intensity of the signal alone) weighted with e.g. efficiency and other parameters for e.g. cross section computation in different data sets simultaneously.

Is there any thinkable way, or is this impossible with RooFit?

Best regards and thanks,
Klaus

[jalopezg]

Hi @canigia; I am sure @moneta can help you with your RooFit question.

Cheers,
J.

[moneta]

Hi,
I don’t think this is easy possible, because if I have well understood, this would require defining the pdf in a complex observables space, and add support for adding using complex coefficients.
A possible solution is for you to code your resulting model pdf and then use that one for fitting in RooFit. You might able to use some of the RooFit components or modify them to adapt for your complex use case

Lorenzo

[canigia]

Hi Lorenzo,

thanks for you feedback!

As alternative I would then try to implement the simultaneous fit myself by using a TMinuit object. Unfortunately I don’t find a more comprehensive tutorial on the web; the Ifit.C example is a bit sparse, i.e. no explanation is given about the meaning of the minuit interface functions and the settings. It doesn’t even explain how to retrieve the fitted parameters or errors after the fit.

Do you (or anybody) know of a more detailed example?

Best,
Klaus

[moneta]

Hi,

The difficult will be to implement your model. You can still use a customised user defined PDF in RooFIt, otherwise you can write yourself the full simultaneous likelihood and minimise using Minuit.
I would not use the TMinuit class, but the ROOT::Fit::Fitter class or the ROOT::Math::Minimiser.
Examples exist in the tutorials/fit directory, like fitCircle.C for the Fitter and NumericalMinimization.C for the Minimiser class.

If you have issues using these classes, please let me know

Best

Lorenzo

[canigia]

Here is my try based on Ifit.C to fit n Gaussians above a parabola background, where mean and sigma of the Gaussians should be constraint to be the same. The printed results, however, look like complete nonsense.

#include "TMinuit.h"

const int nhistmax = 5;
int nhist=0;
TH1F *hfit[nhistmax];

//______________________________________________________________________________
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  double m = par[0];
  double G = par[1];
  
  double chi2 = 0;
  
  // compute chi2 for nhist histograms fitted with Gaus + pol2, Gaus mean and sigma constraint
  for (int ih=0; ih<nhist; ++ih) {
    for (int ib=1; ib<=hfit[ih]->GetNbinsX(); ++ib) {
      int pidx = 2+4*ih;
      
      double x   = hfit[ih]->GetBinCenter(ib);  // x position
      double hy  = hfit[ih]->GetBinContent(ib); // histogram entries
      double hdy = hfit[ih]->GetBinError(ib);   // histogram entries error
      double fy  = par[pidx]*TMath::Gaus(x, m, G, 1) + par[pidx+1] + par[pidx+2]*x + par[pidx+3]*x*x;  // function value
      
      chi2 += hdy>0 ? (hy - fy)/hdy : 0.;
    }
  }
  
  // return fitness
  f = chi2;
}
//______________________________________________________________________________
void minuitfit(int nh = 2, double m = 5, double G = 0.3, int N=10000)
{
  TCanvas *c1 = new TCanvas("c1","c1",20,20,nh*350,350);
  c1->Divide(nh,1,0.0001,0.0001);

  nhist = nh;  
  TF1 *f1 = new TF1("f1","gausn(0) + pol2(3)",0,10);
  
  int npar = 2 + nh*4;
  
  for (int i=0;i<nhist;++i) {
    hfit[i] = new TH1F(Form("h%d",i), Form("hist %d",i), 50+i*50, 0, 10);
    f1->SetParameters(20+i*5, m, G, 20+i*3, 1+i,0.3-0.2*i*i);
    
    hfit[i]->FillRandom("f1",N);
    c1->cd(i+1);
    hfit[i]->Draw("e");
  }
    
  TMinuit *gMinuit = new TMinuit(npar);  //initialize TMinuit with a maximum of 5 params
  gMinuit->SetFCN(fcn);
  Double_t arglist[10];
  Int_t ierflg = 0;
  arglist[0] = 1;
  gMinuit->mnexcm("SET ERR", arglist ,1, ierflg);
  
  // Set starting values and step sizes for parameters
  double vstart[100], step[100];
  
  vstart[0] = m;
  vstart[1] = G;
  for (int i=0; i<nh;++i) {
    vstart[2+i*4]   = 20+i*5;
    vstart[2+i*4+1] = 20+i*3;
    vstart[2+i*4+2] = 1+i;
    vstart[2+i*4+3] = 0.3-0.2*i*i;
  }
  for (int i=0;i<npar; ++i) { 
    step[i] = 0.01;
    gMinuit->mnparm(i, Form("p%02d", i), vstart[i], step[i], 0, 0, ierflg);
  }

  // Print results
  Double_t amin,edm,errdef;
  Int_t nvpar,nparx,icstat;

  // Now ready for minimization step
  arglist[0] = 10000;
  arglist[1] = 0.01;
  gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
  gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
  gMinuit->mnprin(3,amin);
}

[canigia]

Hi Lorenzo,

thanks a lot! I got my simple example now working (there was a bug in the chi^2 computation → the square was missing) both with TMinuit and ROOT::Fit::Fitter.

Best regards,
Klaus