Index: trunk/SophyaExt/MinuitAdapt/Makefile
===================================================================
--- trunk/SophyaExt/MinuitAdapt/Makefile	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/Makefile	(revision 2404)
@@ -0,0 +1,48 @@
+MODULECXXREPNAME :=  MinuitAdapt
+MODULEDECCXXFLAGS := -msg_quiet
+include ../Mgr/Makefile.h
+
+CFLAGS   := $(CFLAGS) -Df2cFortran
+CPPFLAGS := $(CPPFLAGS) -Df2cFortran
+CXXFLAGS := $(CXXFLAGS) -Df2cFortran
+ 
+all: $(LIB)libMinuitAdapt.a
+
+clean: cleantest
+	rm -f $(OBJ)minuitadapt.o $(OBJ)minuitadaptF.o
+	rm -f $(LIB)libMinuitAdapt.a
+
+cleantest:
+	rm -f testminuit testminuit.o
+
+# Doit etre fait automatiquement par mkmflien
+install:
+	cp localmnpout.h minuitprot.h minuitadapt.h $(INC)/.
+
+#----------------------------------------------------------------------
+$(LIB)libMinuitAdapt.a : $(OBJ)minuitadapt.o $(OBJ)minuitadaptF.o
+	$(ARCXX) $(ARCXXFLAGS) $@ $($(ARARGS))
+ifeq ($(CXX),cxx)
+	$(CXX) $? $(CPPFLAGS) $(CXXFLAGS) -o $(OBJ)xx.x -Hf
+endif 
+ 
+$(OBJ)minuitadaptF.o: minuitadaptF.f
+$(OBJ)minuitadapt.o: minuitadapt.cc minuitadapt.h minuitprot.h
+
+#----------------------------------------------------------------------
+ifdef NOSHLIB
+LIBS = -L$(LIB) -lFitsIOServer -lSamba -lSkyMap -lSkyT -lNTools -lHiStats -lTArray -lSysTools -lm
+ifeq ($(CXX),cxx)
+LIBS := $(LIBS) -lertcxx
+endif
+else
+LIBS = -L$(SLB) -lsophya -lextsophya  -lm
+endif
+LIBMINUIT := -L$(LIB) -lMinuitAdapt -L$(EXTLIBPATH) -lminuit
+LDLIBS := $(LIBS) $(LDLIBS) $(LIBMINUIT) $(LIBFORT)
+
+test: testminuit
+
+testminuit: testminuit.o
+	$(LINK.cc) $^ -o $@ $(LOADLIBES) $(LDLIBS)
+$(OBJ)testminuit.o: testminuit.cc
Index: trunk/SophyaExt/MinuitAdapt/README
===================================================================
--- trunk/SophyaExt/MinuitAdapt/README	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/README	(revision 2404)
@@ -0,0 +1,6 @@
+make clean
+make
+make install
+make test
+./testminuit
+make cleantest
Index: trunk/SophyaExt/MinuitAdapt/localmnpout.h
===================================================================
--- trunk/SophyaExt/MinuitAdapt/localmnpout.h	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/localmnpout.h	(revision 2404)
@@ -0,0 +1,17 @@
+/*  Link to minuit (CERNLIB)  */
+/*        CMV - 94            */
+
+#ifndef LOCALMNPOUT_H_SEEN
+#define LOCALMNPOUT_H_SEEN
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void localmnpout_(int_4 *,char *,double *,double *,double *,double *,int_4 *);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
Index: trunk/SophyaExt/MinuitAdapt/minuitadapt.cc
===================================================================
--- trunk/SophyaExt/MinuitAdapt/minuitadapt.cc	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/minuitadapt.cc	(revision 2404)
@@ -0,0 +1,540 @@
+#include "machdefs.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream.h>
+#include <string.h>
+#include <string>
+
+#include "pexceptions.h"
+#include "minuitadapt.h"
+#include "localmnpout.h"
+
+#define AVEC_CFORTRAN
+
+#ifdef AVEC_CFORTRAN
+#include "Cern/cfortran.h"
+#include "Cern/minuit.h"
+#else
+#include "minuitprot.h"
+#endif
+
+/*!
+  \class SOPHYA::Minuit Cern Fitter adapter
+  \ingroup MinuitAdapt
+*/
+
+MinuitAdapt::MinuitAdapt(
+  void (*myfcn)(int_4 *,double *,double *,double *,int_4 *,double (*f)(double *)),
+  double (*myfutils)(double *x)
+  )
+  : mNPar(0), iErFlg(0), fcn(myfcn), futils(myfutils)
+{
+ int_4 i1=5,i2=6,i3=0;
+#ifdef AVEC_CFORTRAN
+ MNINIT(i1,i2,i3);
+#else
+ mninit_(&i1,&i2,&i3);
+#endif
+}
+
+MinuitAdapt::~MinuitAdapt(void)
+{
+}
+
+void MinuitAdapt::SetTitle(char* title)
+{
+ int_4 nc=strlen(title);
+ if(nc>0) {
+#ifdef AVEC_CFORTRAN
+   MNSETI(title);
+#else
+   mnseti_(title,nc);
+#endif
+ }
+}
+
+void MinuitAdapt::Clear(void)
+{
+ char* str="CLEAR";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+ mNPar = 0;
+}
+
+void MinuitAdapt::Return(void)
+{
+ char* str="RETURN";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::PrintLevel(int_4 n)
+{
+ char* str="SET PRINTOUT";
+ double arglis[1]={(double)n}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetWidthPage(int_4 n)
+{
+ char* str="SET WIDTHPAGE";
+ double arglis[1]={(double)n}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetWarnings(bool w)
+{
+ char* str="SET WARNINGS";
+ if(!w) str="SET NOWARNINGS";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetErrorDef(double err)
+// 1=xi2, 0.5=likelyhood
+{
+ char* str="SET ERRORDEF";
+ double arglis[1]={err}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetEpsMachine(double eps)
+{
+ char* str="SET EPSMACHINE";
+ double arglis[1]={eps}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetStrategy(int_4 n)
+{
+ char* str="SET STRATEGY";
+ double arglis[1]={(double)n}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetGradient(int_4 n)
+// n<0 : gradient from finite differences computed by Minuit
+// n=0 : gradient from user AND finite diff.. then compare
+// n>0 : force gradient from user
+{
+ char* str="SET NOGRADIENT";
+ double arglis[1]={1.}; int_4 narg=0;
+ if(n==0) str="SET GRADIENT";
+ else if(n>0) {str="SET GRADIENT"; narg=1;}
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetRandom(int_4 seed)
+// Seed must be integer between 10000 and 900000000
+{
+ char* str="SET RANDOMGENERATOR";
+ double arglis[1]={(double)seed}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::DefineParameter(int_4 n,char* name,double init,double step
+                                 ,double vmin,double vmax)
+// nom parametre, numero parametre, initial value, initial step size, parameter bounds
+// - ATTENTION numero de parametre [1,...]
+{
+ char str1[16]; sprintf(str1,"P%d",n);
+ char* str=str1;
+ if(name!=NULL) str=name;
+#ifdef AVEC_CFORTRAN
+ MNPARM(n,str,init,step,vmin,vmax,iErFlg);
+#else
+ int_4 nc=strlen(str);
+ mnparm_(&n,str,&init,&step,&vmin,&vmax,&iErFlg,nc);
+#endif
+ mNPar++;
+}
+
+void MinuitAdapt::SetParameter(int_4 n,double val)
+{
+ char* str="SET PARAMETER";
+ double arglis[2]={(double)n,val}; int_4 narg=2;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetLimits(int_4 n,double val1,double val2)
+{
+ char* str="SET LIMITS";
+ double arglis[3]={(double)n,val1,val2}; int_4 narg=3;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::SetFix(int_4 n)
+{
+ char* str="FIX";
+ double arglis[1]={(double)n}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Release(int_4 n)
+{
+ char* str="RELEASE";
+ double arglis[1]={(double)n}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Restore(bool only_last_fixed)
+{
+ char* str="RESTORE";
+ double arglis[1]={0.}; int_4 narg=0;
+ if(only_last_fixed) {arglis[0]=1.; narg=1;}
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::DoFit(char *method,int_4 maxcall,double tol)
+{
+ char* str=method;
+ double arglis[2]={(double)maxcall,tol}; int_4 narg=2;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Improve(int_4 maxcall)
+{
+ char* str="IMPROVE";
+ double arglis[1]={(double)maxcall}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Hesse(int_4 maxcall)
+{
+ char* str="HESSE";
+ double arglis[1]={(double)maxcall}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Minos(int_4 maxcall)
+{
+ char* str="MINOS";
+ double arglis[1]={(double)maxcall}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Seek(int_4 maxcall,double stdev)
+{
+ char* str="SEEK";
+ double arglis[2]={(double)maxcall,stdev}; int_4 narg=2;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::Call(int_4 iflag)
+// See meaning in fcn exemple
+{
+ char* str="CALL";
+ double arglis[1]={(double)iflag}; int_4 narg=1;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::DrawContour(int_4 n1,int_4 n2,int_4 npts)
+{
+ char* str="MNCONTOUR";
+ double arglis[3]={(double)n1,(double)n2,(double)npts}; int_4 narg=3;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+int_4 MinuitAdapt::GetContour(int_4 n1,int_4 n2,int_4 npts
+                             ,TVector<r_8>& xpt,TVector<r_8>& ypt)
+{
+ if(npts<5) npts=5;
+ int_4 ncontok = 0;
+ double* xcont = new double[npts];
+ double* ycont = new double[npts];
+#ifdef AVEC_CFORTRAN
+ MNCONT(fcn,n1,n2,npts,xcont[0],ycont[0],ncontok,futils);
+#else
+ mncont_(fcn,&n1,&n2,&npts,xcont,ycont,&ncontok,futils);
+#endif
+ if(ncontok<1) return 0;
+ xpt.ReSize(ncontok); ypt.ReSize(ncontok);
+ for(int_4 i=0;i<ncontok;i++) {xpt(i)=xcont[i]; ypt(i)=ycont[i];}
+ delete [] xcont; delete [] ycont;
+ return ncontok;
+}
+
+void MinuitAdapt::Scan(int_4 n,double from,double to,int_4 npts)
+{
+ char* str="SCAN";
+ double arglis[4]={(double)n,(double)npts,from,to}; int_4 narg=4;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowParameter(void)
+{
+ char* str="SHOW PARAMETER";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowRandom(void)
+{
+ char* str="SHOW RANDOM";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowFcnValue(void)
+{
+ char* str="SHOW FCNVALUE";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowCovariance(void)
+{
+ char* str="SHOW COVARIANCE";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowCorrelations(void)
+{
+ char* str="SHOW CORRELATIONS";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::ShowEigenValues(void)
+{
+ char* str="SHOW EIGENVALUES";
+ double arglis[1]={0.}; int_4 narg=0;
+#ifdef AVEC_CFORTRAN
+ MNEXCM(fcn,str,arglis,narg,iErFlg,futils);
+#else
+ int_4 nc=strlen(str);
+ mnexcm_(fcn,str,arglis,&narg,&iErFlg,futils,nc);
+#endif
+}
+
+void MinuitAdapt::GetParameter(int_4 n,string& name,double& val,double& err
+                              ,double& bound1,double& bound2,int_4& ivarbl)
+{
+ double p,ep,b1,b2; int_4 i; char str[256]="";
+ localmnpout_(&n,str,&p,&ep,&b1,&b2,&i);
+ val=p; err=ep; bound1=b1; bound2=b2; ivarbl=i;
+ // Protection contre les retours intempestifs du Fortran !
+ if(strlen(str)>0)
+   for(unsigned int i=0;i<strlen(str);i++)
+     if(str[i]==' ') {str[i]='\0'; break;}
+ name = str;
+}
+
+void MinuitAdapt::GetErrors(int_4 n,double& eplus,double& eminus
+                           ,double& eparab,double& globcc)
+{
+ double ep,em,epa,g;
+#ifdef AVEC_CFORTRAN
+ MNERRS(n,ep,em,epa,g);
+#else
+ mnerrs_(&n,&ep,&em,&epa,&g);
+#endif
+ eplus=ep; eminus=em; eparab=epa; globcc=g;
+}
+
+TMatrix<r_8> MinuitAdapt::GetErrorsMatrix(void)
+{
+ TMatrix<r_8> mat;
+ int_4 npar=NPar();
+ if(npar<=0) return mat;
+ double* mcov = new double[npar*npar];
+#ifdef AVEC_CFORTRAN
+ MNEMAT(mcov[0],npar);
+#else
+ mnemat_(mcov,&npar);
+#endif
+ mat.ReSize(npar,npar);
+ for(int i=0;i<npar;i++)
+  for(int j=0;j<npar;j++) mat(i,j) = mcov[i*npar+j];
+ delete [] mcov;
+ return mat;
+}
+
+////////////////////////////////////////////////////////////
+/*
+void fcn(int_4 *npar,double *grad,double *fval,double *xval
+        ,int_4 *iflag,double futils(double *))
+{
+  //cout<<"iflag="<<*iflag<<endl;
+  IFLAG[0]++;
+  if(*iflag>0 && *iflag<10) IFLAG[*iflag]++;
+
+  // Read input,init,... data values
+  //  if(*iflag==1) {...}
+
+  // Instruct Minuit to redefine the problem
+  // and forget about previously best fitted values.
+  //  if(*iflag==5) {...}
+
+  // Always compute Chi2 or Likelyhood (here iflag==4)
+ *fval=0.;
+ for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
+   double f = Z[i][j]-Gauss2D(X[i],Y[j],xval);
+   *fval += f*f/(EZ[i][j]*EZ[i][j]);
+ }
+
+ // Compute (optionnal) the first derivative of Chi2 / parameters
+ if(*iflag==2) {
+   // Return gradient of chi2 (if SET GRA called)
+   // C'est DChi2/DPi = -2*sum{(Yi-F(Xi))/EYi^2 * dF/dPi(Xi)}
+   double dpar[NPAR];
+   for(int j=0;j<NPAR;j++) grad[j]=0.;
+   for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
+     double f=-2.*(Z[i][j]-Gauss2D(X[i],Y[j],xval))/(EZ[i][j]*EZ[i][j]);
+     dGauss2D(X[i],Y[j],xval,dpar);
+     for(int k=0;k<NPAR;k++) grad[k]+= f*dpar[k];
+   }
+ }
+
+ // Called at the end of the fit (on the Miniut RETURN)
+ if(*iflag==3) {
+   cout<<"Call fcn iflag="<<*iflag<<endl;
+   for(int k=0;k<NPAR;k++) cout<<" P"<<k+1<<"="<<xval[k];
+   cout<<endl;
+   cout<<"Number of fcn calls="<<IFLAG[0]<<endl;
+   for(int k=1;k<10;k++)
+     cout<<"  iflag="<<k<<" number of calls="<<IFLAG[k]<<endl;
+ }
+}
+*/
Index: trunk/SophyaExt/MinuitAdapt/minuitadapt.h
===================================================================
--- trunk/SophyaExt/MinuitAdapt/minuitadapt.h	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/minuitadapt.h	(revision 2404)
@@ -0,0 +1,79 @@
+// This may look like C code, but it is really -*- C++ -*-
+#ifndef MINUITADAPT_SEEN
+#define MINUITADAPT_SEEN
+
+#include "machdefs.h"
+#include "tvector.h"
+#include "tmatrix.h"
+
+namespace SOPHYA {
+
+class MinuitAdapt {
+public:
+  MinuitAdapt(void (*myfcn)(int_4 *,double *,double *,double *,int_4 *,double (*f)(double *))
+             ,double (*myfutils)(double *)=NULL);
+  virtual ~MinuitAdapt(void);
+
+  void SetTitle(char* title);
+  void Clear(void);
+  void Return(void);
+
+  void PrintLevel(int_4 n=1);
+  void SetWidthPage(int_4 n=120);
+  void SetWarnings(bool w=true);
+  void SetErrorDef(double err=1.);
+  void SetEpsMachine(double eps=1.e-13);
+  void SetStrategy(int_4 n=1);
+  void SetGradient(int_4 n=-1);
+  void SetRandom(int_4 seed);
+
+  void DefineParameter(int_4 n,char* name,double init,double step
+                      ,double vmin=0.,double vmax=0.);
+  void SetParameter(int_4 n,double val);
+  void SetLimits(int_4 n,double val1=0.,double val2=0.);
+  void SetFix(int_4 n);
+  void Release(int_4 n);
+  void Restore(bool only_last_fixed=false);
+
+  void DoFit(char *method="MIGRAD",int_4 maxcall=9999999,double tol=0.1);
+  inline void Migrad(int_4 maxcall=9999999,double tol=0.1)
+       {DoFit("MIGRAD",maxcall,tol);}
+  inline void Minimize(int_4 maxcall=-1,double tol=0.1)
+       {DoFit("MINIMIZE",maxcall,tol);}
+  inline void Simplex(int_4 maxcall=9999999,double tol=0.1)
+       {DoFit("SIMPLEX",maxcall,tol);}
+  void Improve(int_4 maxcall=9999999);
+  void Hesse(int_4 maxcall=9999999);
+  void Minos(int_4 maxcall=9999999);
+  void Seek(int_4 maxcall=9999999,double stdev=3.);
+
+  void Call(int_4 iflag);
+  void DrawContour(int_4 n1,int_4 n2,int_4 npts=25);
+  int_4 GetContour(int_4 n1,int_4 n2,int_4 npts
+                  ,TVector<r_8>& xpt,TVector<r_8>& ypt);
+  void Scan(int_4 n,double from,double to,int_4 npts=25);
+
+  void ShowParameter(void);
+  void ShowRandom(void);
+  void ShowFcnValue(void);
+  void ShowCovariance(void);
+  void ShowCorrelations(void);
+  void ShowEigenValues(void);
+
+  inline int_4 NPar(void) {return mNPar;}
+  inline int_4 ErrFlag(void) {return iErFlg;}
+  void GetParameter(int_4 n,string& name,double& val,double& err
+                   ,double& bound1,double& bound2,int_4& ivarbl);
+  void GetErrors(int_4 n,double& eplus,double& eminus
+                ,double& eparab,double& globcc);
+  TMatrix<r_8> GetErrorsMatrix(void);
+
+protected:
+  int_4 mNPar,iErFlg;
+  void (*fcn)(int_4 *,double *,double *,double *,int_4 *,double (*f)(double *));
+  double (*futils)(double *);
+};
+
+} // Fin du namespace
+
+#endif
Index: trunk/SophyaExt/MinuitAdapt/minuitadaptF.f
===================================================================
--- trunk/SophyaExt/MinuitAdapt/minuitadaptF.f	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/minuitadaptF.f	(revision 2404)
@@ -0,0 +1,12 @@
+*=======================================================
+      subroutine localmnpout(i,c,p,ep,b1,b2,ivar)
+c On surecrit mnpout_:
+c car le C/C++ n'arrive pas a recuperer correctement
+c  la chaine de characteres "c" (dans certains cas core dump)
+      implicit none
+      integer i,ivar
+      real*8 p,ep,b1,b2
+      character *(256) c
+      call mnpout(i,c,p,ep,b1,b2,ivar)
+      return
+      end
Index: trunk/SophyaExt/MinuitAdapt/minuitprot.h
===================================================================
--- trunk/SophyaExt/MinuitAdapt/minuitprot.h	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/minuitprot.h	(revision 2404)
@@ -0,0 +1,32 @@
+/*  Link to minuit (CERNLIB)  */
+/*        CMV - 94            */
+
+#ifndef MINUIT_H_SEEN
+#define MINUIT_H_SEEN
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void mninit_(int_4 *,int_4 *,int_4 *);
+void mnseti_(char *,int_4);
+void mnparm_(int_4 *,char *,double *,double *,double *,double *,int_4 *,int_4);
+void mnexcm_(void fcn (int_4 *,double *,double *,double *,int_4 *,double futil(double *))
+            ,char *,double *,int_4 *,int_4 *,double futil(double *),int_4);
+void mnpout_(int_4 *,char *,double *,double *,double *,double *,int_4 *);
+void mnerrs_(int_4 *,double *,double *,double *,double *);
+void mnemat_(double *,int_4 *);
+void mncont_(void fcn (int_4 *,double *,double *,double *,int_4 *,double futil(double *))
+            ,int_4 *,int_4 *,int_4 *,double *,double *,int_4 *
+            ,double futil(double *));
+
+/*
+void fcn(int_4 *,double *,double *,double *,int_4 *,double futil(double *));
+double futil(double *);
+*/
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
Index: trunk/SophyaExt/MinuitAdapt/testminuit.cc
===================================================================
--- trunk/SophyaExt/MinuitAdapt/testminuit.cc	(revision 2404)
+++ trunk/SophyaExt/MinuitAdapt/testminuit.cc	(revision 2404)
@@ -0,0 +1,309 @@
+/*
+ Fit d'une gaussienne par une methode de chi2
+*/
+#include <machdefs.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <iostream.h>
+#include <math.h>
+#include <string.h>
+
+#include "nbrandom.h"
+#include "minuitadapt.h"
+
+// Ne pas changer NPAR
+#define NPAR 7
+
+const int DIMX = 25;
+const int DIMY = 25;
+const double HAUT=10000.;
+const double X0=0., Y0=0.;
+const double SX=5., SY=5., RHO=0.05;
+const double FOND=100.;
+//                         vol    x0    y0    sx    sy   rho  fond
+const bool FIX[NPAR]={false,false,false,false,false,false,false};
+const double ERR=5.;
+const double ERRMIN=ERR*sqrt(FOND);
+const double nSX=7., nSY=5.;
+const bool USERGRAD=false;
+const bool DOIMPROVE=true;
+const bool DOMINOS=true;
+const bool DOCONT=true;
+const bool DOSCAN=true;
+#define VARAND (drandpm1())
+//#define VARAND 1.
+#define MAXCALL 99999
+
+void fcn(int_4 *,double *,double *,double *,int_4 *,double futils(double *));
+double futils(double *x) {return 0.;}
+double Gauss2D(double x,double y,double *param);
+double dGauss2D(double x,double y,double *param,double *dparam);
+
+double X[DIMX], Y[DIMY], Z[DIMX][DIMY], EZ[DIMX][DIMY];
+int_4 IFLAG[10]={0,0,0,0,0,0,0,0,0,0};
+
+/*==========================================================================*/
+int main(int nargv, char *argv[])
+{
+ /* initialisation */
+ if(DIMX*DIMY<=NPAR) exit(-1);
+
+ int nran=0;
+ if(nargv>1) nran=atoi(argv[1]);
+ for(int i=0;i<nran;i++) drand01();
+
+ /* compute volume */
+ double vol=HAUT*(2*M_PI*SX*SY)/sqrt(1.-RHO*RHO);
+ cout<<"haut="<<HAUT<<" vol="<<vol<<endl;
+ cout<<"Fond="<<FOND<<endl;
+
+ /* remplissage de la densite de probabilite a fiter */
+ double par[NPAR],epar[NPAR];
+ printf("Donnees a fiter %d :\n",DIMX*DIMY);
+   par[0]=vol; par[1]=X0; par[2]=Y0;
+   par[3]=SX; par[4]=SY; par[5]=RHO;
+   par[6]=FOND;
+ for(int i=0;i<DIMX;i++) {
+   X[i] = X0 +nSX*SX*(2.*(i+0.5)/DIMX-1.);
+   for(int j=0;j<DIMY;j++) {
+     Y[j] = Y0 +nSY*SY*(2.*(j+0.5)/DIMY-1.);
+     double f = Gauss2D(X[i],Y[j],par);
+     Z[i][j] = f;
+     double ef = (fabs(f)>1.) ? ERR*sqrt(fabs(f)): ERR;
+     EZ[i][j] = (ef>ERRMIN)? ef: ERRMIN;
+     Z[i][j] += EZ[i][j]*NorRand();
+   }
+ }
+
+/*********************** minuit minimisation ***********************/
+ printf("\n\n");
+ MinuitAdapt MMM(fcn,futils);
+
+ /* initialise */
+ MMM.SetTitle("Minuit fit Gaussienne 2D+Fond");
+ MMM.Clear();
+ MMM.SetRandom(1000000);
+
+ /* set print and warning level, precision etc... (-1,0,1,2,3) */
+ MMM.PrintLevel(1);
+ MMM.SetWidthPage(120);
+ MMM.SetWarnings(true);
+ MMM.SetErrorDef(1.);
+ MMM.SetEpsMachine(1.e-13);
+ MMM.SetStrategy(1);
+ if(USERGRAD) MMM.SetGradient(1); else MMM.SetGradient(-1);
+
+ /* set parameters */
+ MMM.DefineParameter(1,"Vol",vol,fabs(vol)/50.);
+ MMM.DefineParameter(2,"X0",X0,SX/5.);
+ MMM.DefineParameter(3,"Y0",Y0,SY/5.);
+ MMM.DefineParameter(4,"Sx",SX,SX/5.,0.01*SX,10.*SX);
+ MMM.DefineParameter(5,"Sy",SY,SY/5.,0.01*SY,10.*SY);
+ MMM.DefineParameter(6,"Rho",RHO,0.0001,-1.,1.);
+ double x=(FOND!=0.)? fabs(FOND)/10.: 0.01;
+ MMM.DefineParameter(7,"Fond",FOND,x);
+
+ /* set parameters */
+   x = (FIX[0])? vol: vol+VARAND*vol/5.;
+ MMM.SetParameter(1,x);
+   x = (FIX[1])? X0: X0+VARAND*SX;
+ MMM.SetParameter(2,x);
+   x=(FIX[2])? Y0: Y0+VARAND*SY;
+ MMM.SetParameter(3,x);
+   x=(FIX[3])? SX: SX+VARAND*SX/2.;
+ MMM.SetParameter(4,x);
+   x=(FIX[4])? SY: SY+VARAND*SY/2.;
+ MMM.SetParameter(5,x);
+   x=(FIX[5])? RHO: 0.;
+ MMM.SetParameter(6,x);
+   x=(FIX[6])? FOND: FOND+VARAND*FOND/3.;
+ MMM.SetParameter(7,x);
+
+ /* fix parameters */
+ for(int i=0;i<NPAR;i++) if(FIX[i]) MMM.SetFix(i+1);
+ MMM.ShowParameter();
+
+ /* minimize */
+ MMM.Migrad(MAXCALL,0.01);
+ //MMM.Minimize(MAXCALL,0.01);
+ //MMM.Simplex(MAXCALL,0.01);
+ //MMM.Seek(MAXCALL,5.);
+ MMM.ShowFcnValue();
+
+ if(DOIMPROVE) MMM.Improve(MAXCALL);
+ MMM.ShowFcnValue();
+
+ if(DOMINOS) MMM.Minos(MAXCALL);
+ MMM.ShowFcnValue();
+
+ /* get parameters and errors */
+ cout<<endl;
+ for(int i=0;i<NPAR;i++) {
+   string dum; int_4 ivarbl;
+   double b1,b2,eparab,eplus,eminus,globcc;
+   MMM.GetParameter(i+1,dum,par[i],epar[i],b1,b2,ivarbl);
+   MMM.GetErrors(i+1,eplus,eminus,eparab,globcc);
+   printf("> parameter %d \"%s\" = %g %g (%g,%g) int var=%d\n"
+         ,i+1,dum.c_str(),par[i],epar[i],b1,b2,ivarbl);
+   printf("            e+=%g e-=%g eparab=%g globcc=%g\n"
+         ,eplus,eminus,eparab,globcc);
+ }
+ fflush(stdout);
+ cout<<"haut(sig)="<<par[0]*sqrt(1-RHO*RHO)/(2*M_PI*SX*SY)
+     <<"  haut(sig_fit)="<<par[0]*sqrt(1-par[5]*par[5])/(2*M_PI*par[3]*par[4])<<endl;
+
+ TMatrix<r_8> emat = MMM.GetErrorsMatrix();
+ cout<<"GetErrorsMatrix:"<<emat<<endl;
+
+ /* get covariance matrix */
+ MMM.ShowCovariance();
+ MMM.ShowCorrelations();
+ MMM.ShowEigenValues();
+
+ /* contour plot */
+ if(DOCONT) {
+   for(int i=1;i<NPAR;i++) {
+     if(FIX[i]) continue;
+     for(int j=i+1;j<=NPAR;j++) {
+       if(FIX[j]) continue;
+       //MMM.DrawContour(i,j,20);
+       TVector<r_8> xcont,ycont;
+       int_4 ncontok = MMM.GetContour(i,j,20,xcont,ycont);
+       cout<<"Contour "<<i<<" "<<j<<" ncontok="<<ncontok<<endl;
+       if(ncontok<1) continue;
+       for(int k=0;k<ncontok;k++) cout<<" ("<<xcont[k]<<","<<ycont[k]<<")";
+       cout<<endl;
+     }
+   }
+ }
+
+ /* scan parameters */
+ if(DOSCAN) {
+   MMM.ShowFcnValue();
+   for(int i=0;i<NPAR;i++) {
+     if(FIX[i]) continue;
+     MMM.Scan(i+1,par[i]-2.*epar[i],par[i]+2.*epar[i],20);
+   }
+   MMM.ShowFcnValue();
+ }
+
+ /* end */
+ MMM.ShowRandom();
+ MMM.Return();
+
+ printf("\n\n");
+/*******************************************************************/
+
+ /* fin de minimisation, print ! */
+ double xi2=0.;
+ for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
+   double f = Z[i][j]-Gauss2D(X[i],Y[j],par);
+   xi2 += f*f/(EZ[i][j]*EZ[i][j]);
+ }
+ cout<<"1: "<<vol<<"\tfit="<<par[0]<<"\tefit="<<epar[0]<<"\td="<<par[0]-vol;
+   if(epar[0]>0.) cout<<"\td/e="<<(par[0]-vol)/epar[0]; cout<<endl;
+ cout<<"2: "<<X0<<"\tfit="<<par[1]<<"\tefit="<<epar[1]<<"\td="<<par[1]-X0;
+   if(epar[1]>0.) cout<<"\td/e="<<(par[1]-X0)/epar[1]; cout<<endl;
+ cout<<"3: "<<Y0<<"\tfit="<<par[2]<<"\tefit="<<epar[2]<<"\td="<<par[2]-Y0;
+   if(epar[2]>0.) cout<<"\td/e="<<(par[2]-Y0)/epar[2]; cout<<endl;
+ cout<<"4: "<<SX<<"\tfit="<<par[3]<<"\tefit="<<epar[3]<<"\td="<<par[3]-SX;
+   if(epar[3]>0.) cout<<"\td/e="<<(par[3]-SX)/epar[3]; cout<<endl;
+ cout<<"5: "<<SY<<"\tfit="<<par[4]<<"\tefit="<<epar[4]<<"\td="<<par[4]-SY;
+   if(epar[4]>0.) cout<<"\td/e="<<(par[4]-SY)/epar[4]; cout<<endl;
+ cout<<"6: "<<RHO<<"\tfit="<<par[5]<<"\tefit="<<epar[5]<<"\td="<<par[5]-RHO;
+   if(epar[5]>0.) cout<<"\td/e="<<(par[5]-RHO)/epar[5]; cout<<endl;
+ cout<<"7: "<<FOND<<"\tfit="<<par[6]<<"\tefit="<<epar[6]<<"\td="<<par[6]-FOND;
+   if(epar[6]>0.) cout<<"\td/e="<<(par[6]-FOND)/epar[6]; cout<<endl;
+ cout<<"Xi2="<<xi2<<"\txi2/n="<<xi2/(DIMX*DIMY-NPAR)<<endl;
+
+ exit(0);
+}
+
+void fcn(int_4 *npar,double *grad,double *fval,double *xval
+        ,int_4 *iflag,double futils(double *))
+{
+  //cout<<"iflag="<<*iflag<<endl;
+  IFLAG[0]++;
+  if(*iflag>0 && *iflag<10) IFLAG[*iflag]++;
+
+  // Read input,init,... data values
+  //  if(*iflag==1) {...}
+
+  // Instruct Minuit to redefine the problem
+  // and forget about previously best fitted values.
+  //  if(*iflag==5) {...}
+
+  // Always compute Chi2 or Likelyhood (here iflag==4)
+ *fval=0.;
+ for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
+   double f = Z[i][j]-Gauss2D(X[i],Y[j],xval);
+   *fval += f*f/(EZ[i][j]*EZ[i][j]);
+ }
+
+ // Compute (optionnal) the first derivative of Chi2 / parameters
+ if(*iflag==2) {
+   // Return gradient of chi2 (if SET GRA called)
+   // C'est DChi2/DPi = -2*sum{(Yi-F(Xi))/EYi^2 * dF/dPi(Xi)}
+   double dpar[NPAR];
+   for(int j=0;j<NPAR;j++) grad[j]=0.;
+   for(int i=0;i<DIMX;i++) for(int j=0;j<DIMY;j++) {
+     double f=-2.*(Z[i][j]-Gauss2D(X[i],Y[j],xval))/(EZ[i][j]*EZ[i][j]);
+     dGauss2D(X[i],Y[j],xval,dpar);
+     for(int k=0;k<NPAR;k++) grad[k]+= f*dpar[k];
+   }
+ }
+
+ // Called at the end of the fit (on the Minuit RETURN)
+ if(*iflag==3) {
+   cout<<"Call fcn iflag="<<*iflag<<" npar="<<*npar<<endl;
+   for(int k=0;k<NPAR;k++) cout<<" P"<<k+1<<"="<<xval[k];
+   cout<<endl;
+   cout<<"Number of fcn calls="<<IFLAG[0]<<endl;
+   for(int k=1;k<10;k++)
+     cout<<"  iflag="<<k<<" number of calls="<<IFLAG[k]<<endl;
+ }
+}
+
+double Gauss2D(double x,double y,double *param)
+// xc = (x-p1)/p3 ; yc = (y-p2)/p4
+// f(x) = p0*(1-rho^2)/(2*Pi*p3*p4) * exp(-0.5*(xc^2+yc^2-2*p5*xc*yc))
+{
+ double N = sqrt(1.-param[5]*param[5])/(2*M_PI*param[3]*param[4]);
+ double X = (x-param[1])/param[3];
+ double Y = (y-param[2])/param[4];
+ double z2 = (X*X + Y*Y - 2.*param[5]*X*Y)/2.;
+ if(z2<100.) z2=exp(-z2); else z2=0.;
+ return param[0]*N*z2 + param[6];
+}
+
+double dGauss2D(double x,double y,double *param,double *dparam)
+{
+ double unmr2 = 1.-param[5]*param[5];
+ double N = sqrt(unmr2)/(2*M_PI*param[3]*param[4]);
+ double X = (x-param[1])/param[3];
+ double Y = (y-param[2])/param[4];
+
+ double XmrY = X-param[5]*Y;
+ double YmrX = Y-param[5]*X;
+ double z2 = (X*(XmrY-param[5]*Y)+Y*Y)/2.;
+
+ /* g(x,y) */
+ double PSF = 0.;
+ if(z2<100.) PSF = N * exp(-z2);
+ /* dg(x,y)/d(Vol) */
+ dparam[0] = PSF;
+ /* dg(x,y)/d(x0) */
+ dparam[1] =  param[0]* PSF* XmrY/param[3];
+ /* dg(x,y)/d(y0) */
+ dparam[2] =  param[0]* PSF* YmrX/param[4];
+ /* dg(x,y)/d(sx)*/
+ dparam[3] =  param[0]* PSF* (X*XmrY-1.)/param[3];
+ /* dg(x,y)/d(sy) */
+ dparam[4] =  param[0]* PSF* (Y*YmrX-1.)/param[4];
+ /* dg(x,y)/d(rho) */
+ dparam[5] =  param[0]* PSF* (X*Y-2.*param[5]/unmr2);
+ /* dg(x,y)/d(Fond) */
+ dparam[6] = 1.;
+
+ return param[0] * PSF + param[6];
+}
