Problem with minimization (use ROOT::Math::Minimizer)

ROOT
1

Sorry for a lot of code.
I use ROOT’s minimizer for minimizing the following function:

double GKinFitter::Chi2PlusConstraints( const double* unitedArray ){
  
  bool inBarrel;
  chi2 = 0;
 
  double e, theta, phi, rho, z, rhoInitial, zInitial, z0;
  double sigmaE, sigmaP, sigmaR, sigmaZ;
 
  //depending on detection photon either in barrel or edge
  //this variable could have different meaning
  double diffCur;
 
  double dE = constraintE;    //difference between required and measured E
  double dPx = constraintPx;  //... Px 
  double dPy = constraintPy;  //... Py 
  double dPz = constraintPz;  //... Pz 
 
  z0 = unitedArray[ nPars * nPhotons ];//zComInPoint
 
  for( int iPh = 0; iPh < nPhotons; iPh++ ){
    inBarrel = false;
 
    e = unitedArray[ iPh * nPars ];
    sigmaE = photons[iPh]->GetMeasurement(0)->GetError();
 
    phi = unitedArray[ iPh * nPars + 2 ];
    sigmaP = photons[iPh]->GetMeasurement(2)->GetError(); 
 
    rhoInitial = photons[iPh]->GetMeasurement(3)->GetValue();
    sigmaR = photons[iPh]->GetMeasurement(3)->GetError();
    rho = unitedArray[ iPh * nPars + 3 ];
    //check if photon was detected in barrel
    if( rhoInitial > 37.8 ) inBarrel = true;
 
    zInitial = photons[iPh]->GetMeasurement(4)->GetValue();
    sigmaZ = photons[iPh]->GetMeasurement(4)->GetError(); 
    z = unitedArray[ iPh * nPars + 4 ];
 
    if( inBarrel ){
      theta = ArctgPos( rhoInitial / (z + z0) );
      diffCur = pow( (z - zInitial) / sigmaZ, 2 );      
    }
    else{
      theta = ArctgPos( rho / (zInitial + z0) );
      diffCur = pow( (rho - rhoInitial) / sigmaR , 2 );
    }
 
    //Build expressions to minimize (chi2 + constraints)
    chi2 += pow( (e - photons[iPh]->GetMeasurement(0)->GetValue() ) / sigmaE, 2 )\
          + pow( (phi - photons[iPh]->GetMeasurement(2)->GetValue() ) / sigmaP, 2 )\
          + diffCur;
 
    dE -= e;  
    dPx -= e*sin(theta)*cos(phi); 
    dPy -= e*sin(theta)*sin(phi); 
    dPz -= e*cos(theta);  
  }
 
  return chi2 + 5*dE*dE + 5*dPx*dPx + 5*dPy*dPy + 5*dPz*dPz;
 
}

NOTE: This function is a member-function of a class that has one another member function inside which minimizer used:

int GKinFitter::PerformMinimization( const char* minName = "Minuit2",
                                     const char* algoName = "Combined")
{
  //creates the ROOT::Minimizator,
  //sets the minimization parameters ( Tolerance etc.),
  
  int ind;                                //Current index
  double s;                               //Current step
  double ll, ul;                          //Current Limits
  char parName[15];                       //Current Name
 
  double* unitedArray = GatherParameters();
 
  ROOT::Math::Minimizer* minimizer = 
        ROOT::Math::Factory::CreateMinimizer( minName, algoName );
 
  minimizer->SetMaxFunctionCalls( mfcs ); 
  minimizer->SetMaxIterations( mis );
  minimizer->SetTolerance( t ); 
  //Is it correct? 
  ROOT::Math::Functor functor( this, &GKinFitter::Chi2PlusConstraints, ndf );
 
  for ( int iPh = 0; iPh < nPhotons; iPh++ ){ 
    for( int jPar = 0; jPar < nPars; jPar++ ){
 
      ind = iPh*nPars + jPar;
      s = photons[iPh]->GetMeasurement(jPar)->GetError() / sds;
      ll = photons[iPh]->GetMeasurement(jPar)->GetLowerLimit();
      ul = photons[iPh]->GetMeasurement(jPar)->GetUpperLimit();
      sprintf( parName, "freePar%d", ind );
      
      if( jPar != 1 ){  
        minimizer->SetLimitedVariable( ind, parName, unitedArray[ind], s, ll, ul );
      }
      else{
        minimizer->SetFixedVariable( ind, parName, unitedArray[ind] );                          
      }
 
    }
  }
 
  //specify minimization's settings for zComInPoint
  ind = nPars * nPhotons;
  s = zComInPoint.GetError() / sds; 
  ll = zComInPoint.GetLowerLimit(); 
  ul =  zComInPoint.GetUpperLimit(); 
  sprintf(parName,"zComInPoint");
 
  minimizer->SetLimitedVariable( ind, parName, unitedArray[ind], s, ll, ul);
 
  minimizer->SetFunction( functor );
  minimizer->Minimize();  

  delete [] unitedArray;
  delete [] fittedParams;
  minimizer->Clear();
 
  return 0;
 
}

After run this minimization I got the following:

Info in <Minuit2>: Minuit2Minimizer::Minimize : Minimization did NOT converge, Hesse is not valid

And I added in minimized function the following line (in order to monitor changes with variables during minimization (wrong order):

cout << e << " " << phi << " " << rho << " " << z << " " << theta << endl;

And that what I got (numbers are how many times function was called):

2
-2.36 -28.2677 0 40 1.77205 2.186
3
-2.36 -28.2677 0 40 1.77205 2.186
4
-2.36 -28.2677 0 40 1.77205 2.186
5
-2.36 -28.2677 0 40 1.77205 2.186
6
-2.36 -28.2677 0 40 1.77205 2.186
7
-2.36 -28.2677 0 40 1.77205 2.186
8
-2.36 -28.2677 0 40 1.77205 2.186
9
-2.36 -28.2677 0 40 1.77205 2.186
10
-2.36 -28.2677 0 40 1.77205 2.186
11
-2.36 -28.2677 0 40 1.77205 2.186
12
-2.36 -28.2677 0 40 1.77205 2.186
13
-2.36 -28.2677 0 40 1.77205 2.186
14
-2.36 -28.2677 0 40 1.77205 2.186
15
-2.36 -28.2677 0 40 1.77205 2.186
16
-2.36 -28.2677 0 40 1.77205 2.186
17
-2.36 -28.2677 0 40 1.77205 2.186

Why do they stay unchanged?

2

Hi,

The code should be correct. The only think when you create the Functor object pdf, you are passing pdf as the dimension of the functor. Is pdf equals to the total number of parameters (nPhotons * nPars), including the fixed ones ?

The error you are getting is because something is wrong in the function or in the parameter definitions. Probably the derivatives with respects some of the parameters are zero.

I would need the full print out to understand the problem and eventually your full code to run the program

Best Regards

Lorenzo

3

I built more simple function to minimize:

double GKinFitter::Chi2PlusConstraints2( const double* unitedArray ){
                              
  chi2 = 0;
          
  double e, theta, phi;
  double sigmaE, sigmaT, sigmaP;
          
  double dE = constraintE;    //difference between required and measured E
  double dPx = constraintPx;  //... Px 
  double dPy = constraintPy;  //... Py 
  double dPz = constraintPz;  //... Pz 
          
  for( int iPh = 0; iPh < nPhotons; iPh++ ){
          
    e = unitedArray[ iPh * nPars ];
    sigmaE = photons[iPh]->GetMeasurement(0)->GetError();
          
    theta = unitedArray[ iPh * nPars + 1 ];
    sigmaT = photons[iPh]->GetMeasurement(1)->GetError();
          
    phi = unitedArray[ iPh * nPars + 2 ];
    sigmaP = photons[iPh]->GetMeasurement(2)->GetError(); 
          
    //Build expressions to minimize (chi2 + constraints)
    chi2 += pow( (e - photons[iPh]->GetMeasurement(0)->GetValue() ) / sigmaE, 2 )\
          + pow( (theta - photons[iPh]->GetMeasurement(1)->GetValue() ) / sigmaT, 2 )\
          + pow( (phi - photons[iPh]->GetMeasurement(2)->GetValue() ) / sigmaP, 2 );
          
  
    dE -= e;  
    dPx -= e*sin(theta)*cos(phi); 
    dPy -= e*sin(theta)*sin(phi); 
    dPz -= e*cos(theta);      
          
  }
  
  return chi2 + 5*dE*dE + 5*dPx*dPx + 5*dPy*dPy + 5*dPz*dPz;
  
}

Which obviously has minimum (equal to zero). But it is not reached after minimization:

Parameters before minimization: 
0) 1.78205
1) 0.94
2) 0.75
3) 1.76205
4) 2.183
5) -2.36
Function ( x_initial ) = 0.0497828

Parameters after minimization: 
0) 1.76765
1) 0.940252
2) 0.751825
3) 1.75337
4) 2.18311
5) -2.36183
Function ( x_min ) = 0.0385754

Parameters at which it has a minimum should be the following:

Minimum at:
0) 1.73205
1) 0.9553
2) 0.7853
3) 1.73205
4) 2.1863
5) -2.3562

Why the minimum is not reached?

4

I suspect it is connected with tolerance, steps, maxfunccalls and so on. But I am not sure.
By the way, is there a way to print all settings of a minimizer i.e. set variables, its steps, limits; tolerance, … ?

5

Hi,

The values you pass as input are used. You can use a more verbose option, SetPrintLevel(3) and you will see better the output.
But if it does not reach the minimum could be also due to the function itself (maybe a local minimum ).

Lorenzo

6

I set SetPrintLevel(3) and now I think I see ( previous problem):

Info in <Minuit2>: MnHesse: 2nd derivative zero for Parameter  : name = freePar4

I think it goes from this line in objective function:

if( inBarrel ){  
  theta = ArctgPos( rhoInitial / (z + z0) );
  diffCur = pow( (z - zInitial) / sigmaZ, 2 );
}                
else{            
  theta = ArctgPos( rho / (zInitial + z0) );
  diffCur = pow( (rho - rhoInitial) / sigmaR , 2 );
}

if photon has been detected in barrel calorimeter ( in edge calorimeter ) then, indeed, I do not want to variate rho=freePar3 ( z=freePar4 ). My question is
Is it incorrect to do so? I mean should I fix each time this variable? And one more question:
If the variable is fixed but I change it manually in objective function. Is it OK?

7

Hi,

If there is no dependence on that parameter, then you can just declared as fixed, by doing
minimizer->SetFixedVariable( …)
if it is too complicated to remove it from the function

Lorenzo

8

It should be in a function. See why:
Detected photons have some measured parameters:
…, rho, z, …
And depending on where current photon has been detected either rho or z must variates. But each photon has both of this parameters.

9

But I don’t understand. You are computing the chi2 for all the photons (barrel + end-cap) then the function should depends on both rho and z.
If you have minimizing the chi2 for only barrel photons or only end-cap photons you should be able to fix the parameter outside the function.

Lorenzo

10

OK. I got it. If there is no photons in (say) in end-cap then I fix rho. I will try now.

11

It helps, It finishes without Hesse error. But still it is not the minimum. I think that it stops due to some other reason than finding a minimum. By the way, what does +/- mean in the end?:

freePar0	  = 1.87486	 +/-  0.168762	(limited)
freePar1	  = 0.94	 (fixed)
freePar2	  = 0.752884	 +/-  0.0953108	(limited)
freePar3	  = 19.9997	 +/-  2.23909	(limited)
freePar4	  = 14.6067	 (fixed)
freePar5	  = 1.80342	 +/-  0.168763	(limited)
freePar6	  = 2.183	 (fixed)
freePar7	  = -2.36288	 +/-  0.0952927	(limited)
freePar8	  = 19.9997	 +/-  2.09634	(limited)
freePar9	  = -14.0441	 (fixed)
zComInPoint	  = -0.281204	 +/-  0.400662	(limited)

Does it should be read as “The minimum somewhere around this values”?

12

Hi,

How are you so sure is not the minimum ? Have you tried to scan the parameters around the found minimum value.

The uncertainty in the parameters (+/- values) are computed from the inverse of the Hessian. They are the statistical uncertainty in the parameter in the asymptotic assumption that the parameter distributions are normal.

Lorenzo

13

It is obviously not THE minimum because THE minimum is ZERO. I will try to scan the parameters around found.
I am interesting in the following:
You said that

but how to estimate this if I provided only one example of parameter. I.e. sample for each parameter consists of ONE measurement. What sample does minimizer use to estimate sigma?

14

Sorry, Lorenzo. It is not obvious that it is not THE minimum. I forgot the chi2 contribution.
But question about estimation is still valid.

15

Hi,

You have several observable measurements (your data sample), which contribute to the least-square (chi2) sum. From the least-square you find the best estimate of the parameters, but you can also estimate the parameter uncertainties. It is the same when you have a sample of data points, and you can estimate the mean (average of data points value) but you can also estimate the error on the sample mean as the standard deviation / sqrt(N) or sqrt(N-1) if you want the unbiased estimate.

You can find more in any statistics books, for example in Chapter 2 of the book, Data Analysis in High Energy Physics: A Practical Guide to Statistical Methods, or also in the Statistics Chapter of the PDG,
http://pdg.lbl.gov/2016/reviews/rpp2016-rev-statistics.pdf

Lorenzo

16

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.