Home | News | Documentation | Download

Minimize using Minuit and measured values from a Tree

Hi,

I am working with a Tree that contains several properties of my particles, pT, energies, etc. creating a class from it, and working with the .C file.

I am trying to perform a chi2 minimization for each event of my Tree, where I need to minimize a function that looks like this:

chi2 = ( (e1-k1)/sig1 )^2 + ( (e2-k2)/sig2 )^2 + ((M - 2k1k2*cost)/gam )^2,

where e1,e2, cost are measured valued that I can obtain from the Tree.
M, gam, sig1, sig2 are constant values defined inside my ::Loop(), and k1,k2 are the parameters I want to fit to make the chi2 as close to zero as possible.

The questions are:

  • First, I am not sure where should I be defining the fcn function. It doesn’t let me do it inside the myclass::Loop, so up to now I have been doing it before the myclass:Loop function.

  • And second: I need to use some measured variables from the Tree to enter into the chi2 function and I don’t know how to do it. I am just declaring the parameters to fit i.e. k1, and k2, but don’t find a way to define e1,e2,cost.

Part of the code I have got up to now:

void fcn(int &npar, double *gin, double &f, double *par, int iflag)
{
   // e1,e2,cost needed from Tree. 
   //M,gam,sig1,sig2 constants defined inside myclass:Loop.
   double k1 = par[0];
   double k2 = par[1];
   double Mk12 = 2 * k1 * k2 * (1-cost);
   f = pow( (e1-k1) / sig1 ,2) + pow( (e2-k2) / sig2 ,2) + pow( (M- Mk12)/gam ,2)
}

void myclass::Loop(){
if (fChain == 0) return;
Long64_t nentries = fChain->GetEntriesFast();
Long64_t nbytes = 0, nb = 0;

for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

      double sig1=15., sig2=15., gam=2., M=100.;
      double e1=jte[iq1]; //measured values from my Tree
      double e2=jte[iq2];
      double cost = jte[iq3]; 

      TMinuit *gMinuit = new TMinuit(2); 
      gMinuit->SetFCN(fcn);
      gMinuit->DefineParameter(0, "k1", e1 , 1., 0., 0.); //initial estimates equal to measured values
      gMinuit->DefineParameter(1, "k2", e2 , 1., 0., 0.);
      gMinuit->Command("MIGRAD");
      double x,y,xerr,yerr; 
      gMinuit->GetParameter(0,x,xerr); 
      gMinuit->GetParameter(1,y,yerr);

      printf("x: %+.7f +- %.7f\n",x,xerr);
      printf("y: %+.7f +- %.7f\n",y,yerr);
      }
 }

And of course when I try to run this I get error like the following:

**error:** **use of undeclared identifier 'e1'**

Thanks for any help.

ROOT Version: 6.22/02
Built for macosx64 on Aug 17 2020, 12:46:52
From tags/v6-22-02@v6-22-02

1 Like

Hello,
How big is your tree ? If the values you are retrieving from the Tree can fit easly in memory, I would first retrieve these data, put in some arrays (e.g. an std::vector) and then later perform the loop where you compute the FCN function.
I would also not use directly TMinuit, but the ROOT::Math::Minimizer interface where you can easly define a Functor class for your FCN and pass to the Minimizer.
See ROOT: tutorials/fit/NumericalMinimization.C File Reference for an example on how to use the class

Lorenzo

Total number of events: 325320

But I will need to use this code for another Trees that may have even larger size.

Thanks for the ROOT::Math::Minimizer. But, again how can I create my function depending on the measured variables? In the example the RosenBrock only has the two parameters to fit as arguments. The two main lines should be these, but I’m not sure how to express this dependence.

double RosenBrock(const double *xx ) //definition of my function to minimize

ROOT::Math::Functor f(&RosenBrock,2); //creating the Functor class from my defined function.

Thanks!

Hi,
The number of entries is not too large and you can put them in memory. This will make the computation of your chi2 function much faster.
Suppose your variables are "e1" and "e2". You can do:

Long64_t nentries = fChain->GetEntriesFast();
fChain->Draw("e1:e2"); 
std::vector<double> e1(fChain->GetV1(), fChain->GetV1()+nentries); 
std::vector<double> e2(fChain->GetV1(), fChain->GetV1()+nentries); 
auto myFCN = [&](const double * par) {
    double chi2 = 0;  
    for (int i = 0; i < nentries; i++) {
       double k1 = par[0];
      double k2 = par[1];
     chi2 += pow( (e1[i]-k1) / sig1 ,2) + pow( (e2[i]-k2) / sig2 ,2);
 } 
 return chi2; 
};
ROOT::Math::Functor f(myFCN,2); // 2 is the number of dimension of the FCN 

and use f later as in the tutorial example to minimise the function.

Lorenzo

1 Like

Hi Lorenzo, thanks a lot!

I ended also trying to just doing it without creating these arrays (std::vector) by doing everything inside my Loop. In part because some of the variables I needed were not exactly a measured value, but a calculated value from measured values (after applying some selections).

Here is how it looks. Just a final doubt for setting some limits at the end.

//Before this part I have defined jte1, jte2, cost, sig1,sig2,gam
 
auto myFCN = [&](const double * par) {
     double chi2 = 0;  
     double k1 = par[0];
     double k2 = par[1];
     double M12 = sqrt(2.*k1*k2*(1.-cost));
     chi2 = (jte1-k1)*(jte1-k1)/sig1/sig1 + (jte2-k2)*(jte2-k2)/sig2/sig2 + (M12-MW)*(M12-MW)/gam ;
     return chi2; 
  };
  ROOT::Math::Functor f(myFCN,2);
     //-------------------------------------------
     //-----------------------------------------------------------
     ROOT::Math::Minimizer* minimum =
        ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
  
     double step[2] = {0.001,0.001};
     // starting point
     double variable[2] = { jte[iq1],jte[iq2]}; //same as measured value

     minimum->SetFunction(f);

     // Set the free variables to be minimized !
     //minimum->SetVariable(0,"jte1",variable[0], step[0]);
     //minimum->SetVariableLowerLimit(0, 0.);

     minimum->SetLowerLimitedVariable(0,"k1",variable[0], step[0],0); //Lower value is energy=0
     minimum->SetLowerLimitedVariable(1,"k2",variable[1], step[1],0);
     
     // do the minimization
     minimum->Minimize();

     const double *xs = minimum->X();
     cout << "Minimum: f(" << xs[0] << "," << xs[1] << "): "
              << minimum->MinValue()  << std::endl;

I tried to set some limits on my fit variables, but it is giving me a warning which I don’t understand.

I tried to set the lower limit of k1 to 0, as I can’t get negative values of energy, but it returns the following:

Warning in <TMinuitMinimizer::SetLowerLimitedVariable>: not implemented - use as upper limit 1.E7 instead of +if

Not sure why it says anything about upper limit (and +inf) if I’m trying to set a lower limit, and I could’t find some implementation of this online. Neither the SetVariableLowerLimit works. So, any help is appreciated.

Thanks!

Try with: "Minuit2"

1 Like

Thanks @Wile_E_Coyote . I solved the warning using Minuit2, but now I have some cases where the minimization does not converges. For that event I get the following:

Minuit2Minimizer: Minimize with max-calls 420 convergence for edm < 0.01 strategy 1
MnSeedGenerator: for initial parameters FCN = 10132.67344111
MnSeedGenerator: Initial state: - FCN = 10132.67344111 Edm = 26793.7 NCalls = 9
VariableMetric: start iterating until Edm is < 2e-05
VariableMetric: Initial state - FCN = 10132.67344111 Edm = 26793.7 NCalls = 9
VariableMetric: Iteration # 0 - FCN = 10132.67344111 Edm = 26793.7 NCalls = 9
Info: VariableMetricBuilder: no improvement in line search
VariableMetric: Iteration # 1 - FCN = 10132.67344111 Edm = 26793.7 NCalls = 20
Info: VariableMetricBuilder: iterations finish without convergence.
Info in VariableMetricBuilder : edm = 107175
Info in requested : edmval = 2e-05
Info in matrix forced pos-def by adding to diagonal : padd = 0.0488126
Info: MnHesse: matrix was forced pos. def.
VariableMetric: After Hessian - FCN = 10132.67344111 Edm = 42744.1 NCalls = 30
VariableMetric: Iteration # 2 - FCN = 10132.67344111 Edm = 42744.1 NCalls = 30
Info: VariableMetricBuilder: Tolerance is not sufficient, continue the minimization
Info in Current Edm is : edm = 42744.1
Info in Required Edm is : edmval = 2e-05
Info: DavidonErrorUpdator: delgam < 0 : first derivatives increasing along search line
VariableMetric: Iteration # 3 - FCN = 9908.358418474 Edm = 45773.4 NCalls = 37
Info: VariableMetricBuilder: no improvement in line search
VariableMetric: Iteration # 4 - FCN = 9908.358418474 Edm = 45773.4 NCalls = 48
Info: VariableMetricBuilder: iterations finish without convergence.
Info in VariableMetricBuilder : edm = 1.5399e+06
Info in requested : edmval = 2e-05
Info: VariableMetricBuilder: FunctionMinimum is invalid after second try
Info: MnGlobalCorrelationCoeff: inversion of matrix fails.
Number of iterations 5
----------> Iteration 0
FVAL = 10132.6734411 Edm = 26793.6755668 Nfcn = 9
Error matrix change = 1
Parameters : p0 = 266.617 p1 = 67.4984
----------> Iteration 1
FVAL = 10132.6734411 Edm = 26793.6755668 Nfcn = 20
Error matrix change = 1
Parameters : p0 = 266.617 p1 = 67.4984
----------> Iteration 2
FVAL = 10132.6734411 Edm = 42744.0554093 Nfcn = 30
Error matrix change = 1
Parameters : p0 = 266.617 p1 = 67.4984
----------> Iteration 3
FVAL = 9908.35841847 Edm = 45773.3953658 Nfcn = 37
Error matrix change = 10.8806
Parameters : p0 = 273.237 p1 = 61.0848
----------> Iteration 4
FVAL = 9908.35841847 Edm = 45773.3953658 Nfcn = 37
Error matrix change = 10.8806
Parameters : p0 = 273.237 p1 = 61.0848
Info in Minuit2Minimizer::Minimize : Minimization did NOT converge, Edm is above max
Minuit2Minimizer : Invalid Minimum - status = 3
FVAL = 9908.36
Edm = 45773.4
Nfcn = 37

Here is the piece of relevant code:

These are some values obtained before minimization.
jte[ijb] = 136.945
cost = 0.156321
vtop = 35696.2
jte[iq1] = 266.617
jte[iq2] = 67.4984
ptx= 49.822
pty=6.5479
ptz=-363.926

  double MW=80.379;
  double Mtop=173.3, gammatop = 1.35;
  double sig1=15., sig2=15., gammaW=2.085;

  auto myFCN = [&](const double * par) {
     double chi2 = 0;  
     double k1 = par[0];
     double k2 = par[1];
     double M12 = sqrt(2.*k1*k2*(1.-cost));
     double kb = 2*(k1+k2)*jte[ijb];
     double mtjj = sqrt(M12+kb-vtop); //top mass from Wjets + bjet
     double etoph = k1+k2+jte[ijb] ; // Energy of hadronic Top Reco
     double mtoph = sqrt(etoph*etoph - ptx*ptx - pty*pty - ptz*ptz) ;

     chi2 = (jte[iq1]-k1)*(jte[iq1]-k1)/sig1/sig1 + (jte[iq2]-k2)*(jte[iq2]-k2)/sig2/sig2   
              + (M12-MW)*(M12-MW)/gammaW2 + (mtoph-Mtop)*(mtoph-Mtop)/gammatop/gammatop ;

     return chi2; 
  };
  ROOT::Math::Functor f(myFCN,2);
  ROOT::Math::Minimizer* minimum = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
  minimum->SetTolerance(0.01);
  minimum->SetPrintLevel(3);
  minimum->SetFunction(f);

  double step[2] = {0.01,0.01};
  double variable[2] = { jte[iq1],jte[iq2] }; // starting point

  minimum->SetLowerLimitedVariable(0,"k1",variable[0], step[0],0.); //Lower value is energy=0
  minimum->SetLowerLimitedVariable(1,"k2",variable[1], step[1],0.);
   minimum->Minimize(); // do the minimization

Reproducible example. minimisho.C (1.7 KB)

Thanks!

So, I tried two cases with “partial chi2”:

  • chi2 = (M12-MW)*(M12-MW)/gammaW/gammaW; returns k1 = 123 and k2 = 31
  • chi2 = (mtoph-Mtop)*(mtoph-Mtop)/gammatop/gammatop; returns k1 = 207 and k2 = 62

Note that the returned best fit results differ roughly by a factor of 2 ("k1 * k2" ratio is about 3.4 and “k1 + k2” ratio is about 1.75). The second case returns values not that distant from the initial values (the first case finds a minimum roughly at MW = initial M12 / 2).
To me, this suggests that you might need to review the equations which calculate “M12” and “mtoph” (otherwise, a common minimum will not really be found).
Also, the analytical “chi2” formula in your first post above is a sum of 3 terms, while your newest source code contains 4 terms (so maybe you should review the way you calculate “chi2”).

Thanks.

Why when using “Minuit” in this case it does find the minima? In that case k1=254.087, and k2=15.1852 and a reasonable value for the chi2. I still use the SetLowerLimitedVariable to define my initial params. If I only use SetVariable it again returns error.

Is there a way to use “Minuit” and set the lower limits of my parameters? Or another way to fix it?

Chi2 in my first post was the simplest case. When I was sure the minimizer was working for the simple case I added the extra term at the end.

Thanks!

Well, minimizers in ROOT return a “local minimum”. There is no guarantee that a “global minimum” will be found if the function has many “local minima”.

It is crucial then that the initial values of all fitted parameters are as close as possible to the desired minimum.

Note that “M12” depends only on “k1*k2” (whose initial value is around 18000). You got k1=254 and k2=15 so “k1*k2=3810” and this corresponds to my first test case, which found the same “k1*k2” value (my second test case returned a significantly different “k1*k2=12800”, which is much closer to the initial value, but “mtoph” depends on “k1+k2” only, which should then be used for comparisons).