Get detector info from fitting MC to data

[jingliu]

I need to shift and smear energies obtained from MC so that the resulting distribution matches that of data. I wrote the following code to demonstrate the idea. However, the fitting never converges. How can I solve the problem? Thanks, Jing

TH1F *h0;
double shift0, smear0, scale0, e[3000];

double model(double *x, double *p)
{        
  double shift = p[0];
  double smear = p[1];
  double scale = p[2];
                                                                                        
  if (shift!=shift0 || smear!=smear0 || scale!=scale0) {
    shift0=shift; smear0=smear; scale0=scale; // save old values
    delete h0; h0 = new TH1F("h0","", 50, -5, 15); // empty histogram
    // shift, smear and then scale distribution of e[3000]
    for (int i=0; i<3000; i++) h0->Fill(gRandom->Gaus(e[i]*shift, smear), scale);
  }
  return h0->GetBinContent(int(2.5*(x[0]+5))+1);
}

void test()
{
  // create fake data histogram
  TF1 fgaus("fgaus","exp(-0.5*((x-9)/1.2)**2)",4,15);
  TH1F *h1 = new TH1F("h1", "", 50, -5, 15);
  h1->FillRandom("fgaus", 100);

  // create a set of fake MC samples
  for (int i=0; i<3000; i++) e[i]=gRandom->Gaus(45,1);

  // create model using fake MC samples
  TF1 *f = new TF1("f", model, 4, 15, 3);
  f->SetParNames("shift", "smear", "scale");
  f->SetParameters(0.2,1.2,0.1); 
                                                                                        
  // fit model to data
  gStyle->SetOptFit(1011);
  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2","Simplex");
  ROOT::Math::MinimizerOptions::PrintDefault("Minuit2");

  h1->Fit(f,"rlv");
  h0->SetLineColor(kBlue);
  h0->Draw("same");
}

[couet]

Yes, I see it does not converge. May be the initial parameters are too far from the final values / I guess @moneta can help.

[jingliu]

I’ve tried different initial values, adding limits, different minuit algorithms, tolerance, etc. They help to some extent, the fitting result looks better. But I never managed to get it converge. The fitting is too sensitive to initial values and the fitting range. I am wondering if I did something fundamentally wrong.

[Wile_E_Coyote]

In principle, fitting algorithms expect that the fitting function and its (first two?) derivatives with respect to parameters behave well. That’s not your case. You explicitly introduce significant “random jitter” (“random noise”) in it.

[jingliu]

@Wile_E_Coyote, thanks for the tip. I changed return h0->GetBinContent(int(2.5*(x[0]+5))+1); to return h0->Interpolate(x[0]); to avoid sudden changes in the function. I also smoothed the function 5 times using h0->Smooth(5,"R"). The function looks pretty smooth to me by eye:

But the fitting still cannot converge.

[Wile_E_Coyote]

In your “model” function, try:

  // ...
  double v = 0.;
  for (unsigned long i = 0; i < (sizeof(e) / sizeof(double)); i++)
    v += TMath::Gaus(x[0], e[i] * shift, smear, 1); // 0 or 1
  return scale * v;

[jingliu]

@Wile_E_Coyote, it worked! Thank you!

By the way, I got

Info in Minuit2Minimizer::Minimize : Covar is not pos def
Minuit2Minimizer : Valid minimum - status = 5

The covariance and correlation matrix elements are all zero. I don’t expect the parameters to be correlated. But do you know why all the elements are zero? And what does status 5 mean? Where can I find the explanation of status values?

Thanks, Jing

[Wile_E_Coyote]

See: ROOT::Minuit2::Minuit2Minimizer::Minimize

The “status = 5” is triggered by the “Covar is not pos def” (@moneta “Hesse” and/or “MnSeed” failed and then all elements of the covariance and correlation matrices are zero, too?)

It is also possible that it’s because you use the "Simplex" minimization algorithm (which possibly does not produce the covariance matrix at all). Try, e.g.:

ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2", "Migrad"); // "Migrad" or "Combined"