An example of use of the minimization class TMinuit

//*CMZ :  0.90/09 03/12/96  17.45.55  by  Rene Brun
//*-- Author :    Rene Brun   22/08/95
// ---------------------------------- minexam.C

// Copyright (C) 1994 CodeCERN. All rights reserved.

//______________________________________________________________________________
//*-*-*-*-*-*-*-*-*-*-*-*The Minuit standard test program-*-*-*-*-*-*-*-*-*
//*-*                    ========================                         *
//*-*                                                                     *
//*-*    This program is the translation to C++ of the minexam program    *
//*-*    distributed with the Minuit/Fortran source file.                 *
//*-*         original author Fred James                                  *
//*-*                                                                     *
//*-*       Fit randomly-generated leptonic K0 decays to the              *
//*-*       time distribution expected for interfering K1 and K2,         *
//*-*       with free parameters Re(X), Im(X), DeltaM, and GammaS.        *
//*-*                                                                     *
//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

#include "TROOT.h"
#include "TMinuit.h"
#include "TMath.h"
#include "TStopwatch.h"

#include <stdlib.h>
#include <stdio.h>
#include <iostream.h>

#ifdef __SC__
long G__globalvarpointer; // To make the Symantec linker happy
#endif

int Error;

extern void fcnk0(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);

static ncount = 0;

VoidFuncPtr_t dictfuncs[] = { NULL };

//______________________________________________________________________________
main()
{
   TROOT minexam("TestMinuit","Test of Minuit", dictfuncs);

   static Int_t nprm[5] = { 0   ,    1   ,     4    ,   9     ,  10};
   static Double_t vstrt[5] = {0.  ,    0.  ,    .535  ,   .892   ,  518.3};
   static Double_t stp[5] = {0.1 ,    0.1 ,     0.1  ,     0.   ,   0.};
   static Double_t p0=0;
   static Double_t p1=1;
   static Double_t p3=3;
   static Double_t p5=5;
   Int_t ierflg = 0;
   const Int_t FIVE=5;
   const Int_t THREE=3;
   TStopwatch timer;

   TMinuit *gMinuit = new TMinuit(5);  //initialize TMinuit with a maximum of 5 params
   cout<<"Starting timer"<<endl;
   timer.Start();
   gMinuit->mninit(5,6,7);
   gMinuit->SetFCN(fcnk0);

   gMinuit->mnparm(nprm[0], "Re(X)",    vstrt[0], stp[0], 0,0,ierflg);
   gMinuit->mnparm(nprm[1], "Im(X)",    vstrt[1], stp[1], 0,0,ierflg);
   gMinuit->mnparm(nprm[2], "Delta M",  vstrt[2], stp[2], 0,0,ierflg);
   gMinuit->mnparm(nprm[3], "T Kshort", vstrt[3], stp[3], 0,0,ierflg);
   gMinuit->mnparm(nprm[4], "T Klong",  vstrt[4], stp[4], 0,0,ierflg);
   if (ierflg) {
      Printf(" UNABLE TO DEFINE PARAMETER NO.");
      return ierflg;
   }

//   gMinuit->mnseti("Time Distribution of Leptonic K0 Decays");
//*-*-       Request FCN to read in (or generate random) data (IFLAG=1)
   gMinuit->mnexcm("CALL FCN", &p1 ,1,ierflg);

   gMinuit->mnexcm("FIX", &p5 ,1,ierflg);
   gMinuit->mnexcm("SET PRINT", &p0 ,1,ierflg);
   gMinuit->mnexcm("MIGRAD", &p0 ,0,ierflg);
   gMinuit->mnexcm("MINOS", &p0 ,0,ierflg);
   gMinuit->mnexcm("RELEASE", &p5 ,1,ierflg);
   gMinuit->mnexcm("MIGRAD", &p0 ,0,ierflg);
   gMinuit->mnexcm("MINOS",  &p0 ,0,ierflg);
   gMinuit->mnexcm("CALL FCN", &p3 , 1,ierflg);
   cout<<"Time at the end of job = "<<timer.CpuTime()<<" seconds"<<endl;
}

//______________________________________________________________________________
void fcnk0(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag)
{
   const Int_t MXBIN=50;
   static Double_t thplu[MXBIN],thmin[MXBIN],t[MXBIN];
   static Int_t nbins = 30;
   static Int_t nevtot = 250;
   static Double_t evtp[30] = {
              11.,  9., 13., 13., 17.,  9.,  1.,  7.,  8.,  9.,
               6.,  4.,  6.,  3.,  7.,  4.,  7.,  3.,  8.,  4.,
               6.,  5.,  7.,  2.,  7.,  1.,  4.,  1.,  4.,  5.};
   static Double_t evtm[30] = {
               0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.,  1.,
               0.,  2.,  1.,  4.,  4.,  2.,  4.,  2.,  2.,  0.,
               2.,  3.,  7.,  2.,  3.,  6.,  2.,  4.,  1.,  5.};

   static Int_t i, nevplu, nevmin;
   static Double_t xre, xim, dm, gams, gaml,gamls,xr1,xr2,em,ep;
   static Double_t sthplu, sthmin, ehalf, th, sterm, sevtp, sevtm;
   static Double_t thplui, thmini, chisq, ti, thp, thm, evp, evm, chi1, chi2;

   xre   = x[0];
   xim   = x[1];
   dm    = x[4];
   gams  = 1/x[9];
   gaml  = 1/x[10];
   gamls = 0.5*(gaml+gams);
   if (iflag != 1)  goto L300;
//                        generate random data
   sthplu = sthmin = 0;
   for (i=1;i<=nbins; i++) {
      t[i-1] = 0.1*Double_t(i);
      ti = t[i-1];
      ehalf = TMath::Exp(-ti*gamls);
      xr1 = 1-xre;
      xr2 = 1+xre;
      th =      (xr1*xr1 + xim*xim) * TMath::Exp(-ti*gaml);
      th = th + (xr2*xr2 + xim*xim) * TMath::Exp(-ti*gams);
      th = th -               4*xim*TMath::Sin(dm*ti) * ehalf;
      sterm = 2*(1-xre*xre-xim*xim)*TMath::Cos(dm*ti) * ehalf;
      thplu[i-1] = th + sterm;
      thmin[i-1] = th - sterm;
      sthplu += thplu[i-1];
      sthmin += thmin[i-1];
   }
   nevplu = Int_t(Double_t(nevtot)*(sthplu/(sthplu+sthmin)));
   nevmin = Int_t(Double_t(nevtot)*(sthmin/(sthplu+sthmin)));
   cout<<"  LEPTONIC K ZERO DECAYS"<<endl;
   cout<<Form(" PLUS, MINUS, TOTAL=%5d %5d %5d",nevplu,nevmin,nevtot)<<endl;
   cout<<"0      TIME       THEOR+      EXPTL+      THEOR-      EXPTL-"<<endl;
   sevtp = sevtm = 0;
   for (i=1;i<=nbins; i++) {
      thplu[i-1] = thplu[i-1]*Double_t(nevplu) / sthplu;
      thmin[i-1] = thmin[i-1]*Double_t(nevmin) / sthmin;
      thplui = thplu[i-1];
      sevtp += evtp[i-1];
      thmini = thmin[i-1];
      sevtm += evtm[i-1];
      if (iflag != 4) {
         cout<<Form("%12.4f%12.4f%12.4f%12.4f%12.4f",t[i-1],thplu[i-1],evtp[i-1],thmin[i-1],evtm[i-1])<<endl;
      }
   }
   cout<<Form(" DATA EVTS PLUS, MINUS= %10.2f%10.2f", sevtp,sevtm)<<endl;

//                      calculate chisquared
L300:
   chisq = sthplu = sthmin = 0;
   for (i=1;i<=nbins; i++) {
      ti = t[i-1];
      ehalf = TMath::Exp(-ti*gamls);
      xr1 = 1-xre;
      xr2 = 1+xre;
      th =      (xr1*xr1 + xim*xim) * TMath::Exp(-ti*gaml);
      th = th + (xr2*xr2 + xim*xim) * TMath::Exp(-ti*gams);
      th = th -               4*xim*TMath::Sin(dm*ti) * ehalf;
      sterm = 2*(1-xre*xre-xim*xim)*TMath::Cos(dm*ti) * ehalf;
      thplu[i-1] = th + sterm;
      thmin[i-1] = th - sterm;
      sthplu += thplu[i-1];
      sthmin += thmin[i-1];
  }
   thp = thm = evp = evm = 0;
   if (iflag != 4) cout<<"          POSITIVE LEPTONS                    ,NEGATIVE LEPTONS"<<endl;

   if (iflag != 4) {
      cout<<"      TIME    THEOR    EXPTL    chisq         TIME    THEOR    EXPTL    chisq"<<endl;
   }
   for (i=1;i<=nbins; i++) {
      thplu[i-1] = thplu[i-1]*sevtp / sthplu;
      thmin[i-1] = thmin[i-1]*sevtm / sthmin;
      thp += thplu[i-1];
      thm += thmin[i-1];
      evp += evtp[i-1];
      evm += evtm[i-1];
//  Sum over bins until at least four events found
      if (evp > 3)  {
         ep   = evp-thp;
         chi1 = (ep*ep)/evp;
         chisq = chisq + chi1;
         if (iflag != 4) {
            cout<<Form(" %9.3f%9.3f%9.3f%9.3f",t[i-1],thp,evp,chi1)<<endl;
         }
         thp = evp = 0;
      }
      if (evm > 3) {
         em   = evm-thm;
         chi2 = (em*em)/evm;
         chisq = chisq + chi2;
         if (iflag != 4) {
            cout<<Form("                                          %9.3f%9.3f%9.3f%9.3f"
                             ,t[i-1],thm,evm,chi2)<<endl;
         }
         thm = evm = 0;
      }
  }
      f = chisq;
      ncount++;
//      cout<<ncount<< " F= "<<f<<" xre="<<xre<<" xim="<<xim<<" dm="<<dm<<" gams="<<gams<<" gaml="<<gaml<<endl;
}