Changeset 392 in Sophya for trunk

Timestamp:

Aug 13, 1999, 11:57:20 PM (26 years ago)

Author:

ercodmgr

Message:

Fitter avec init auto cmv 13/8/99

Location:

trunk/SophyaPI/PIext

Files:

• piafitting.cc (modified) (18 diffs)
• piafitting.h (modified) (2 diffs)

trunk/SophyaPI/PIext/piafitting.cc

@@ -6 +6 @@
 //clic sur FIT 2sd fois -> core dump!!!
 //
-//
 #include <stdio.h>
 #include <stdlib.h>
 #include <ctype.h>
-
 #include <iostream.h>
-
 #include <typeinfo>
 
-#include "piafitting.h"
 #include "strutil.h"
-
+#include "nbtri.h"
 #include "generalfit.h"
-
 #include "fct1dfit.h"
 #include "fct2dfit.h"
-
 #include "matrix.h"
 #include "cvector.h"
 #include "ntuple.h"
 #include "cimage.h"
-
 #include "histos.h"
 #include "histos2.h"
@@ -33 +26 @@
 #include "hisprof.h"
 
+#include "piafitting.h"
 #include "nobjmgr.h"
 #include "pistdimgapp.h"
@@ -179 +173 @@
 , mFit(NULL)
 , mV(NULL), mH(NULL), mM(NULL), mH2(NULL), mIm(NULL), mG(NULL)
-, mFunction(NULL), mFName(""), mFunc(NULL), mFitUserFunc(false)
+, mFunction(NULL), mFName(""), mFunc(NULL)
 , mDlUFunc(NULL), mNameFitFunc(""), mFitFunc(NULL), mFitFuncDer(NULL)
 , mUFNVar(0), mUFNPar(0)
@@ -507 +501 @@
    return;}
 if(mFunction!=NULL) {delete mFunction; mFunction=NULL;}
-if(mFunc!=NULL)     {delete mFunc;     mFunc=NULL;}
-mNPar=0;
-mFitUserFunc = false;
+if(mFunc    !=NULL) {delete mFunc;     mFunc=NULL;}
+mNPar=0; mFName = "";
 
 if(func == mNameFitFunc) { // Fonction utilisateur
   mFunc = new GeneralFunc(mUFNVar,mUFNPar,mFitFunc,mFitFuncDer);
   mNPar = mUFNPar;
-  mFitUserFunc = true;
   cout<<"Fonction utilisateur nvar="<<mUFNVar<<", npar="<<mNPar<<endl;
 } else if(func[0]=='p' && mNVar==1) { //polynome
@@ -556 +548 @@
   cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
   if(integ) {GauRhInt2D* myf = new GauRhInt2D; mFunction = myf;}
-    else    {GauRho2D* myf = new GauRho2D; mFunction = myf;}
+    else    {GauRho2D*   myf = new GauRho2D;   mFunction = myf;}
   mNPar = mFunction->NPar(); mFName = func;
 
@@ -563 +555 @@
   cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
   if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; mFunction = myf;}
-    else    {GdlRho2D*   myf = new GdlRho2D; mFunction = myf;}
+    else    {GdlRho2D*   myf = new GdlRho2D;   mFunction = myf;}
   mNPar = mFunction->NPar(); mFName = func;
 
@@ -570 +562 @@
   cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
   if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; mFunction = myf;}
-    else    {Gdl1Rho2D*   myf = new Gdl1Rho2D; mFunction = myf;}
+    else    {Gdl1Rho2D*   myf = new Gdl1Rho2D;   mFunction = myf;}
   mNPar = mFunction->NPar(); mFName = func;
 
@@ -577 +569 @@
   cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
   if(integ) {MofRhInt2D* myf = new MofRhInt2D; mFunction = myf;}
-    else    {MofRho2D*   myf = new MofRho2D; mFunction = myf;}
+    else    {MofRho2D*   myf = new MofRho2D;   mFunction = myf;}
    mNPar = mFunction->NPar(); mFName = func;
 
@@ -606 +598 @@
 // Initialisation automatique des parametres du fit
 {
-cout<<"CMV: to be done : PIAFitter::AutoIniFit"<<endl;
-if(mFitUserFunc) { // Fonction utilisateur
+if(mFunc) { // Fonction utilisateur
   cout<<"PIAFitter::AutoIniFit : Fonction utilisateur, pas d'init auto"<<endl;
-
 } else if(mFName[0]=='p' && mNVar==1) { //poly 1D
   cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 1D"<<endl;
-
-} else if(mFName[0]=='e' && mNVar==1) { //expo+pol 1D
-  cout<<"PIAFitter::AutoIniFit : Auto Init type expo+pol 1D"<<endl;
+  IniFitP1ou2D();
 
 } else if(mFName[0]=='g' && mNVar==1) { //gauss en haut 1D
   cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_haut+pol 1D"<<endl;
+  IniFitGhP1D();
 
 } else if(mFName[0]=='G' && mNVar==1) { //gauss en vol 1D
   cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+pol 1D"<<endl;
+  IniFitGvP1D();
 
 } else if(mFName[0]=='p' && mNVar==2) { //poly 2D
   cout<<"PIAFitter::AutoIniFit : Auto Init type polynome 2D"<<endl;
+  IniFitP1ou2D();
 
 } else if((mFName[0]=='G'||mFName[0]=='d') && mNVar==2) { //gauss+fond en vol
   cout<<"PIAFitter::AutoIniFit : Auto Init type gauss_vol+fond 2D"<<endl;
+  IniFitGv2D();
 
 } else if(mFName[0]=='D' && mNVar==2) { //DLgauss+fond en vol 2D
   cout<<"PIAFitter::AutoIniFit : Auto Init type DLgauss_vol+fond 2D"<<endl;
+  IniFitGv2D();
   Vector dum(mNPar);
   dum = mPar;  mPar(6)=0.5;   mPar(7)=dum(6);
@@ -639 +632 @@
 } else if(mFName[0]=='M' && mNVar==2) { //Moffat+fond en vol 2D
   cout<<"PIAFitter::AutoIniFit : Auto Init type moffat_vol+fond 2D"<<endl;
+  IniFitGv2D();
   Vector dum(mNPar);
   dum = mPar;  mPar(6)=2.5;   mPar(7)=dum(6);
@@ -715 +709 @@
 // classe GeneraFit
 if(mFit != NULL) {delete mFit; mFit= NULL;}
-if(mFunction) mFit= new GeneralFit(mFunction);
-  else        mFit= new GeneralFit(mFunc);
+if(mFunction)    mFit = new GeneralFit(mFunction);
+  else if(mFunc) mFit = new GeneralFit(mFunc);
+if(!mFit) return(-1002);
+
 // Options et Set de GeneralFit
 mFit->SetDebug(mOpt.lpg);
@@ -767 +763 @@
     if(mOpt.okfun) {Matrix* ob = mM->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
   } else if(mH2) {
-    if(mOpt.okres) {Histo2D* ob = mH2->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
-    if(mOpt.okfun) {Histo2D* ob = mH2->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
+    if(mOpt.okres)
+      {Histo2D* ob = mH2->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
+    if(mOpt.okfun)
+      {Histo2D* ob = mH2->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
   } else if(mIm) {
     if(mOpt.okres) {
@@ -779 +777 @@
     }
   } else if(mG) {
-    if(mOpt.okres) {GeneralFitData* ob = mG->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
-    if(mOpt.okfun) {GeneralFitData* ob = mG->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
+    if(mOpt.okres)
+      {GeneralFitData* ob = mG->FitResidus(*mFit);  if(ob) omg.AddObj(ob,nomres);}
+    if(mOpt.okfun)
+      {GeneralFitData* ob = mG->FitFunction(*mFit); if(ob) omg.AddObj(ob,nomfun);}
   }
 }
@@ -923 +923 @@
   return(-1);
 }
+
 // decodage des arguments de la commande
 string mSp,mSs,mSm,mSM,mSf,mSo;
@@ -938 +939 @@
     else if(c[0]=='o') {mSo += ","; mSo += c+2;}
   }
+
 // Execution des commandes
 DecodeOptions(mSo);
@@ -960 +962 @@
 }
 
+/* --Methode-- */
+void PIAFitter::IniFitP1ou2D(void)
+{
+if(mNPar==0) return;
+mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
+if(!mGData) return;
+double mean,sigma;
+int rc = mGData->GetMeanSigma(0,mean,sigma);
+if(rc<=0) return;
+mPar(0) = mean;
+return;
+}
+
+void PIAFitter::IniFitGhP1D(void)
+{
+if(mNPar==0) return;
+mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
+if(mNVar!=1 || mNPar<3) return;
+double h,v,x0,y0,sx,sy,f;
+int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
+if(rc!=0) return;
+mPar(0)=h;  mStep(0)=h/1000.;
+mPar(1)=x0; mStep(1)=sx/10.;
+mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx; 
+if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=h+f;
+}
+
+void PIAFitter::IniFitGvP1D(void)
+{
+if(mNPar==0) return;
+mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
+if(mNVar!=1 || mNPar<3) return;
+double h,v,x0,y0,sx,sy,f;
+int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
+if(rc!=0) return;
+mPar(0)=v;  mStep(0)=v/1000.;
+mPar(1)=x0; mStep(1)=sx/10.;
+mPar(2)=sx; mStep(2)=sx/10.; mMin(2)=0.; mMax(2)=20.*sx; 
+if(mNPar>3) {mPar(3)=f; mStep(3)=f/1000.;} else mPar(0)=v+mNData*f;
+}
+
+void PIAFitter::IniFitGv2D(void)
+{
+if(mNPar==0) return;
+mPar=0.; mStep=1.; mMin=1.; mMax=-1.; mFix=0.;
+if(mNVar!=2 || mNPar<6) return;
+double h,v,x0,y0,sx,sy,f;
+int rc = IniFitGaus(*mGData,h,v,x0,y0,sx,sy,f);
+if(rc!=0) return;
+mPar(0)=v;  mStep(0)=v/1000.;
+mPar(1)=x0; mStep(1)=sx/10.;
+mPar(2)=y0; mStep(2)=sy/10.;
+mPar(3)=sx; mStep(3)=sx/10.; mMin(3)=0.;  mMax(3)=20.*sx; 
+mPar(4)=sy; mStep(4)=sy/10.; mMin(4)=0.;  mMax(4)=20.*sy; 
+mPar(5)=0.; mStep(5)=0.005;  mMin(5)=-1.; mMax(5)=1.; 
+if(mNPar>6) {mPar(6)=f; mStep(6)=f/1000.;} else mPar(0)=v+mNData*f;
+}
+
+int PIAFitter::IniFitGaus(GeneralFitData& g,double& h,double& v
+    ,double& x0,double& y0,double& sx,double& sy,double& f)
+// Pour initialiser une gaussienne 1D ou 2D
+// return code =0 si OK
+{
+h=v=x0=y0=sx=sy=f=0.;
+if(!&g) return(-1);
+int nvar = g.NVar();
+int ndata = g.NData();
+if(nvar<1 || 2<nvar || ndata<=0) return(-2);
+
+// X et Y ordonnees croissant
+double *X=NULL, *Y=NULL;
+X = new double[ndata];
+if(nvar==2) Y = new double[ndata];
+{for(int i=0;i<ndata;i++) {X[i] = g.X(i); if(Y) Y[i] = g.Y(i);}}
+qsort(X,(size_t) ndata,sizeof(double),qSort_Dble);
+if(Y) qsort(Y,(size_t) ndata,sizeof(double),qSort_Dble);
+
+// determination des zones externes=[0,1/3] et [2/3,1]
+//            et des zones internes=[1/3,2/3]
+int k1=ndata/3,k2=2*ndata/3;
+double xmin,xmax, ymin=-1.,ymax=1.;
+xmin = X[k1]; xmax = X[k2];
+if(Y) {ymin = Y[k1]; ymax = Y[k2];}
+
+// recherche des parametres
+double nb=0, nc=0;
+{for(int i=0;i<ndata;i++) {
+   double val,x,y=0.;
+   x = g.X(i); if(Y) y = g.Y(i); val = g.Val(i);
+   if(xmin<x && x<xmax && ymin<y && y<ymax) { // c'est le centre
+     if(nc==0 || val>h) {h=val; x0=x; y0=y;}
+     v += val; nc++;
+   } else {f += val; nb++;} // c'est le bord
+
+}}
+if(nb==0 || nc==0) {delete [] X; if(Y) delete [] Y; return(-3);}
+
+// Calcul des parametres
+f /= (double) nb; h -= f; v -= f*nc;
+sx = (xmax-xmin)/2.; if(Y) sy = (ymax-ymin)/2.; 
+
+delete [] X; if(Y) delete [] Y;
+return(0);
+}
+
+
 ////////////////////////////////////////////////////////////////////
 ////---------- Classe de fenetre de fit interactive ------------////
@@ -995 +1103 @@
 cpx = spx; cpy = spy;
 string dum = "Object: "+mFitter->mNObj+" , Fun: ";
-if(mFitter->mFitUserFunc) dum+=mFitter->mNameFitFunc; else dum+=mFitter->mFName;
+if(mFitter->mFunc) dum+=mFitter->mNameFitFunc; else dum+=mFitter->mFName;
 lab[0] = new PILabel(this,dum.c_str(),4*bsx,bsy,cpx,cpy);
 lab[0]->SetBinding(PIBK_elastic,PIBK_elastic, PIBK_elastic,PIBK_elastic);
@@ -1351 +1459 @@
 }
 else if(2011<=msg && msg<=2013) {
-  mFitter->mOpt.polcx = msg-2011;
+  mFitter->mOpt.polcy = msg-2011;
   if(lp) cout<<"Centrage Y polcy="<<mFitter->mOpt.polcy
              <<" yc="<<mFitter->mOpt.yc<<endl;

trunk/SophyaPI/PIext/piafitting.h

@@ -69 +69 @@
 void LinFit(void);
 
+// pour initialiser les parametres du fit
+void IniFitP1ou2D(void);
+void IniFitGhP1D(void);
+void IniFitGvP1D(void);
+void IniFitGv2D(void);
+static int IniFitGaus(GeneralFitData& g,double& h,double& v
+       ,double& x0,double& y0,double& sx,double& sy,double& f);
+
 // Graphique
 PIStdImgApp* mApp; // Do not delete
@@ -89 +97 @@
 GeneralFunction* mFunction; string mFName;
 // Gestion de fonction de fit definie par l'utilisateur
-bool mFitUserFunc;
 GeneralFunc* mFunc;
 PDynLinkMgr* mDlUFunc;