Pointwise function in minuit

Hi
I would like to use minuit to find the matrix elements of a standard model calculation, by comparing measured vs simulated histograms. In all the examples I have seen, minuit needs a TF1 or TF2 to compare to the histogram. But I am comparing the output of a monte carlo simulation to the measurement data.

Instead of using a TF2 directly, is there a way of comparing measured data to a histogram weighted by a TF2 and then extracting the weights/parameters? In other words, I cannot define my model as a pure TF2 but must define it as another histogram weighted by a TF2.

Many thanks!

Hi,

You can define the model function you are using to fit also as another histograms. You need just to include in that function, which is based on an histogram, some parameters which represents a continuous variation.
A possible parameter for example is the overall histogram normalisation. So the shape is fixed but you are fitting the histogram normalisation to the observed data.
Another example, at the opposite extreme, is to define a different parameter for each the bin content of the histogram. In this case you will overfit the data, you have tot many parameters and the resulting chi-square will be always zero.
Technically you can implement a TF2 using a user defined function, or a user class (a functor), which can be built from whatever object you need.

Best Regards

Lorenzo

At present my (stripped down) program looks like the below. Is this how I would account for the histogram? If I simply switch fparam to a TH2F throughout, then operations like fparam->SetParameter(i,par[i]) will no longer be defined.

Note: the function “matrixCalc” is very long and not included.

TH2F *hdata;
TH2F *hsim;
TF2  *fparam;

double calcNLL(TH2F* h_meas, TH2F* h_sim, TF2 *f){
  double nll=0;
  TAxis *xaxis = h_meas->GetXaxis();
  TAxis *yaxis = h_meas->GetYaxis();
  for (int i=1; i<=h_meas->GetNbinsX(); i++)for (int j=1; j<=h_meas->GetNbinsY(); j++){
      double X=xaxis->GetBinCenter(i);
      double Y=yaxis->GetBinCenter(j);
      int n=(int)(h_meas->GetBinContent(i,j));
      double mu=(double)h_sim->GetBinContent(i,j);
      mu=mu*f->Eval(X,Y);
      if (mu<1e-10) mu=1e-10;   
      nll += n * TMath::Log(mu) + mu  + TMath::LnGamma(n+1);
    }
  return -2*nll;  
}
void fcn(int &npar, double *deriv, double &f, double par[], int flag){
  for (int i=0; i<npar; i++){
    fparam->SetParameter(i,par[i]);
  }
  f = calcNLL(hdata,hsim,fparam);
} 
main(int argc, char **argv) {
   TH2F *hexp=(TH2F*)expFile->Get("X_Y_meas");
   TF2* myfunc = new TF2("myfunc", matrixCalc, 0, 1, 0, 1, npar);  

   hdata=hexp;
   fparam=myfunc;
  TMinuit minuit(npar);
   minuit.SetFCN(fcn); 
   minuit.Migrad();       
   double outpar[npar], err[npar];
   for (int i=0; i<npar; i++){
     minuit.GetParameter(i,outpar[i],err[i]);
    }
}