Understand MINUIT

Dil · June 16, 2022, 7:40pm

Hello Experts,
I am having some trouble understanding how to use MINUIT class for minimization function.
I have the following function
func= TMath::Power(x,par[0]) * TMath::Power((1-x),par[1])TMath::Power((1+ (yy)/ par[3]),-par[2]);
I was thinking to use LogLikelihood method for minimization(hoping to use MINUIT package for parameter minimization)

I have tree variables called xBj ,QSq…etc.I have divided 2D plot into the nQq=4 and nXB=3.What I need to do is,
need to find the minimization function for each Q2 X xB bin.
So that I defined TF2 for each bin. I am having trouble understanding how I minimize the function using MINUIT for different parameters.
I have attached my code here.

execute the following code

.L xBQ2binning.cpp
Could you please help me to understand how to use MINUIT in this case?

I tried to follow the following example.But I am not sure to continue this.
https://root.cern/doc/master/Ifit_8C.html

Apology if this is a weird question.

Thank You
Dil

xBQ2binning_tree.h (4.0 KB)
xBQ2binning.cpp (4.9 KB)
xBQ2binning_tree.C (1.9 KB)

xBQ2binning.root (1.9 MB)

moneta · June 17, 2022, 2:55pm

Hi,

I think you would like to fit an histogram with the given function defined as a TF2. In this case you don’t need to worry too much about the minimization algorithm, but just do:

TF2 f2 = new TF2("f2","TMath::Power(x,par[0]) * TMath::Power((1-x),par[1])TMath::Power((1+ (yy)/ par[3]),-par[2])");
h2->Fit(f2);

Instead if you need to minimise a function yourself, like a log-likelihood function, then it is different.
In this case first implement a function class and then you can use the Minimizer interface to minimise it.
See the tutorial, NumericalMinimization.C.

You can also use directly the legacy interface of TMinuit, as in the tutorial you have linked above, but I think is more complex.

Lorenzo

Dil · June 22, 2022, 1:59pm

Hello Lorenzo,

As you suggested I tried to follow your example. But I am having following crash.This crash happens when I introduced the following minimization step.

minimum[jQ2Bin][ixBin]->Minimize();

Here is my code,

//
//  xBQ2binning.cpp
//
//
//
//

#include <stdio.h>
#define xBQ2binning_tree_cxx
#include "xBQ2binning_tree.h"
#include "xBQ2binning_tree.C"

#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"

//   In a ROOT session, you can do:
//      root> .L xBQ2binning_tree.C
//      root> xBQ2binning_tree t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//double fcn(double *par);
//void myFCN(int &npar, double *gin, double &logli, double *par, int iflag);
void xBQ2binning(){
    // contained in .h file for this Class
    //TFile *input =new TFile("xBQ2binning_simu.root","read");
    //  Put event sum of logli function here for just one jQ2bin, ixBin value
    int jQ2,ixB;
    int nbin=0;
    double Norm[nQ2][nxB], PhSp[nQ2][nxB];
    double  NBin[nQ2][nxB];
    TRandom3 ran3;
    double xBj,QSq;
    TTree          *fChain;
    
    ROOT::Math::Minimizer* minimum[nQ2][nxB];
    //xBQ2binning_tree tntuple;
    TMinuit *gMinuit[nQ2][nxB];
    int nparMINUIT=4;
    TNtuple *xBQ2binning_tree = new TNtuple("xBQ2binning_tree","Branches","xBj:QSq:ixB:jQ2");
    char fName[64];
    TF2 *fLikeli[nQ2][nxB];
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            sprintf(fName,"fLikeli_Q2_%02d_xB_%02d",jQ2Bin, ixBin);
            fLikeli[jQ2Bin][ixBin]= new TF2(fName,Likelihood,0.0,1.0,0.0,10.0,4);
            // revise to different parameter sets for each Q2 x xB bin
            fLikeli[jQ2Bin][ixBin]->SetParameters(-0.5,2.0,2.0,1.0);
        }
    }
    
    
    double xBMinBin, xBMaxBin;
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            Norm[jQ2Bin][ixBin] = 0.0;
            PhSp[jQ2Bin][ixBin]= 0.0;
            for (int iSimu=0; iSimu<nsimu; iSimu++){
                QSq = QSqBinLow[jQ2Bin]+ (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin])*ran3.Uniform();
                xBMinBin = (xBjMin(QSq)*(nxB-ixBin)+ xBjMax(QSq)*ixBin )/(double)nxB;
                xBMaxBin = (xBjMin(QSq)*(nxB-ixBin-1)+ xBjMax(QSq)*(ixBin+1) )/(double)nxB;
                xBj = xBMinBin+(xBMaxBin-xBMinBin)*ran3.Uniform();
                Norm[jQ2Bin][ixBin] += (fLikeli[jQ2Bin][ixBin]->Eval(xBj,QSq))*(xBMaxBin-xBMinBin);
                PhSp[jQ2Bin][ixBin] += (xBMaxBin-xBMinBin);
            }
            Norm[jQ2Bin][ixBin] *= (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin])/nsimu;
            PhSp[jQ2Bin][ixBin] *= (QSqBinLow[jQ2Bin+1]-QSqBinLow[jQ2Bin]) / nsimu;
            
        }
    }
    //    Put event loop inside bin loop for minimization purposes
    for (int jQ2Bin=0; jQ2Bin<nQ2; jQ2Bin++) {
        for (int ixBin=0; ixBin<nxB; ixBin++) {
            
            
            
            
            
            auto myFCN = [&](const double *par){
                
                double logli=0.0;
                int jQ2Bin = par[4];
                int ixBin  = par[5];
              //  if (fChain == 0) return;
                Long64_t nentries = fChain->GetEntries();
             //   Long64_t nentries = fChain->GetEntriesFast();
                fChain->GetEntriesFast();
                // int nentries= xBQ2binning_tree->GetEntries();
                
                Long64_t nbytes = 0, nb = 0;
                for (Long64_t jentry=0; jentry<nentries;jentry++) {
                    Long64_t ientry = fChain->LoadTree(jentry);  //Long64_t ientry = LoadTree(jentry);
                    
                    if (ientry < 0) break;
                    nb = fChain->GetEntry(jentry);   nbytes += nb;
                    
                    // Select events from bin (ixBin, jQ2Bin)
                    if ((ixB==ixBin)&&(jQ2==jQ2Bin)) {
                        // replace with explicit function
                        logli += log(TMath::Power(xBj,par[0]) * TMath::Power((1-xBj),par[1])*TMath::Power((1+ (QSq*QSq)/ par[3]),-par[2]))-log(Norm[jQ2Bin][ixBin]); //not sure about the normalization is correct here
                        nbin++;
                    }
                }
                logli *= -1.0/nbin;
                //       printf("LogLi[jQ2=%02d,ixB=%02d]= %10.3f \n", jQ2Bin, ixBin, logli);
                
                return logli;
                
            };
            ROOT::Math::Functor f(myFCN,2);
            minimum[jQ2Bin][ixBin] = ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
            minimum[jQ2Bin][ixBin]->SetTolerance(0.01);
            minimum[jQ2Bin][ixBin]->SetPrintLevel(3);
            minimum[jQ2Bin][ixBin]->SetFunction(f);
            
            double step[6] = {0.2 , 0.2 , 0.2 , 0.2, 0.01, 0.01};
            double variable[6] = { -0.5, 2.0 , 2.0 , 1.0, (double)jQ2Bin, (double)ixBin }; // starting point
            
            
            minimum[jQ2Bin][ixBin]->SetVariable(0,"alpha",variable[0], step[0]);
            minimum[jQ2Bin][ixBin]->SetVariable(1,"beta",variable[1], step[1]);
            minimum[jQ2Bin][ixBin]->SetVariable(2,"gamma",variable[2], step[2]);
            minimum[jQ2Bin][ixBin]->SetVariable(3,"lambda",variable[3], step[3]);
            minimum[jQ2Bin][ixBin]->SetVariable(4,"jQ2Bin",variable[4], step[4]);
            minimum[jQ2Bin][ixBin]->SetVariable(5,"ixBin",variable[5], step[5]);
          //  minimum[jQ2Bin][ixBin]->SetFixedVariable(3);
            //  minimum[jQ2Bin][ixBin]->SetFixedVariable(4);
             //minimum[jQ2Bin][ixBin]->SetFixedVariable(5);
            
            
            // do the minimization
            minimum[jQ2Bin][ixBin]->Minimize();
            //   minimum[jQ2Bin][ixBin]->SetLowerLimitedVariable(0,"k1",variable[0], step[0],0.); //Lower value is energy=0
            // minimum[jQ2Bin][ixBin]->SetLowerLimitedVariable(1,"k2",variable[1], step[1],0.);
            //minimum[jQ2Bin][ixBin]->Minimize(); // do the minimization
            
        }
        
        
    }
}

Crash

**********
 **    1 **SET PRINT           2
 **********
 **********
 **    2 **SET NOGRAD
 **********
 PARAMETER DEFINITIONS:
    NO.   NAME         VALUE      STEP SIZE      LIMITS
     1 alpha       -5.00000e-01  2.00000e-01     no limits
     2 beta         2.00000e+00  2.00000e-01     no limits
     3 gamma        2.00000e+00  2.00000e-01     no limits
     4 lambda       1.00000e+00  2.00000e-01     no limits
     5 jQ2Bin       0.00000e+00  1.00000e-02     no limits
     6 ixBin        0.00000e+00  1.00000e-02     no limits
 **********
 **    3 **SET ERR           1
 **********
 **********
 **    4 **SET PRINT           2
 **********
 **********
 **    5 **SET STR           1
 **********
 NOW USING STRATEGY  1: TRY TO BALANCE SPEED AGAINST RELIABILITY
 **********
 **    6 **MIGRAD           0        0.01
 **********
 FIRST CALL TO USER FUNCTION AT NEW START POINT, WITH IFLAG=4.
[/Users/Dil/root-6-12-04-02/lib/libCling.so] cling_runtime_internal_throwIfInvalidPointer (no debug info)
[<unknown binary>] (no debug info)
[<unknown binary>] (no debug info)
[<unknown binary>] (no debug info)
[<unknown binary>] (no debug info)
[/Users/Dil/root-6-12-04-02/lib/libMinuit.so] TMinuitMinimizer::Fcn(int&, double*, double&, double*, int) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuitMinimizer.cxx:253
[/Users/Dilini/root-6-12-04-02/lib/libMinuit.so] TMinuit::Eval(int, double*, double&, double*, int) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuit.cxx:830
[/Users/Dilini/root-6-12-04-02/lib/libMinuit.so] TMinuit::mnamin() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuit.cxx:992
[/Users/Dilini/root-6-12-04-02/lib/libMinuit.so] TMinuit::mnmigr() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuit.cxx:5072
[/Users/Dilini/root-6-12-04-02/lib/libMinuit.so] TMinuit::mnexcm(char const*, double*, int, int&) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuit.cxx:2859
[/Users/Dilini/root-6-12-04-02/lib/libMinuit.so] TMinuitMinimizer::Minimize() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/math/minuit/src/TMinuitMinimizer.cxx:562
[<unknown binary>] (no debug info)
[<unknown binary>] (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] cling::IncrementalExecutor::executeWrapper(llvm::StringRef, cling::Value*) (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] cling::Interpreter::RunFunction(clang::FunctionDecl const*, cling::Value*) (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] cling::Interpreter::EvaluateInternal(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, cling::CompilationOptions, cling::Value*, cling::Transaction**, unsigned long) (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] cling::Interpreter::process(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> > const&, cling::Value*, cling::Transaction**, bool) (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] cling::MetaProcessor::process(llvm::StringRef, cling::Interpreter::CompilationResult&, cling::Value*, bool) (no debug info)
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] HandleInterpreterException(cling::MetaProcessor*, char const*, cling::Interpreter::CompilationResult&, cling::Value*) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/metacling/src/TCling.cxx:2087
[/Users/Dilini/root-6-12-04-02/lib/libCling.so] TCling::ProcessLine(char const*, TInterpreter::EErrorCode*) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/metacling/src/TCling.cxx:2236
[/Users/Dilini/root-6-12-04-02/lib/libRint.so] TRint::ProcessLineNr(char const*, char const*, int*) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/rint/src/TRint.cxx:741
[/Users/Dilini/root-6-12-04-02/lib/libRint.so] TRint::HandleTermInput() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/rint/src/TRint.cxx:603
[/Users/Dilini/root-6-12-04-02/lib/libCore.so] TUnixSystem::CheckDescriptors() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/unix/src/TUnixSystem.cxx:1323
[/Users/Dilini/root-6-12-04-02/lib/libCore.so] TMacOSXSystem::DispatchOneEvent(bool) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/macosx/src/TMacOSXSystem.mm:378
[/Users/Dilini/root-6-12-04-02/lib/libCore.so] TSystem::InnerLoop() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/base/src/TSystem.cxx:412
[/Users/Dilini/root-6-12-04-02/lib/libCore.so] TSystem::Run() /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/base/src/TSystem.cxx:362
[/Users/Dilini/root-6-12-04-02/lib/libCore.so] TApplication::Run(bool) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/base/src/TApplication.cxx:1208
[/Users/Dilini/root-6-12-04-02/lib/libRint.so] TRint::Run(bool) /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/core/rint/src/TRint.cxx:458
[/Users/Dilini/root-6-12-04-02/bin/root.exe] main /Users/Dilini/Documents/codetest/GEMC_Simulation/montecarloPiplusminus/root/main/src/rmain.cxx:32
[/usr/lib/system/libdyld.dylib] start (no debug info)
[<unknown binary>] (no debug info)
Error in <HandleInterpreterException>: Trying to access a pointer that points to an invalid memory address..
Execution of your code was aborted.
In file included from input_line_10:1:
/Users/Dilini/Documents/Maximum_likeli/Maximum_likeli/Testing_minimization.cpp:99:37: warning: invalid memory pointer passed to a callee:
                Long64_t nentries = fChain->GetEntries();
                                    ^~~~~~
root [2]

Attached the files. Could you please help me to figure this out?

xBQ2binning.root (1.9 MB)
Testing_minimization.cpp (6.3 KB)
xBQ2binning_tree.C (1.9 KB)
xBQ2binning_tree.h (4.0 KB)

moneta · June 22, 2022, 2:51pm

It seems you are not assigning the fChain pointer, and it therefore null during the minimization

Lorenzo

Dil · June 22, 2022, 3:03pm

Thank you Lorenzo.But I was wondering why do I need to define fChain again in my main code because I already defined in .h file.
So I just declare it locally(inside cpp file) as follows
TTree *fChain;
Thanks
Dil

moneta · June 22, 2022, 3:30pm

You have defined as a data member of the class and now you are using in a function. These are two separate scopes.
You need to declare and assign it to the TTree you load from the file and you need to be sure it is alive during the minimization time when it is used in your lambda function.

Lorenzo