Changeset 361 in Sophya for trunk


Ignore:
Timestamp:
Aug 6, 1999, 7:16:22 PM (26 years ago)
Author:
ercodmgr
Message:

Deplacement des methodes d'ajustement dans nouvelle classe Reza 6/8/99

Location:
trunk/SophyaPI/PIext
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/SophyaPI/PIext/basexecut.cc

    r357 r361  
    566566  }
    567567
    568 
    569 // Fit 1D sur objets 1D. Egalement Fit 2D sur objets 2D.
    570 else if (kw == "fit") {
    571   if (tokens.size() < 2) {
    572     cout <<"Usage:fit nomobj func \n"
    573          <<" [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]\n";
    574     return(0);
    575   }
    576   string p=""; string s=""; string m=""; string M=""; string O="";
    577   if (tokens.size()>2)
    578     for(int ip=2;ip<tokens.size();ip++) {
    579       if(tokens[ip].length()<=2) continue;
    580       const char *c = tokens[ip].c_str();
    581       if(c[1]!=':') continue;
    582       if(c[0]=='p')      p=c+2;
    583       else if(c[0]=='s') s=c+2;
    584       else if(c[0]=='m') m=c+2;
    585       else if(c[0]=='M') M=c+2;
    586       else if(c[0]=='o') {O += ","; O += c+2;}
    587     }
    588   srvo->Fit12D(tokens[0],tokens[1],p,s,m,M,O);
    589 }
    590568
    591569
     
    918896mpiac->RegisterCommand(kw, usage, this, "Expr. Plotting");
    919897
    920 kw = "fit";
    921 usage = "Fitting function to DataObjects (Histo, Histo2D, Vector, ...)";
    922 usage += "\n Usage: fit nomobj func [Options]";
    923 usage += "\n [p:p1,...,pn s:s1,...,sn m:m1,...,mn M:M1,...,Mn o:... o:...]";
    924 mpiac->RegisterCommand(kw, usage, this, "Expr. Plotting");
    925 
    926898}
    927899
  • trunk/SophyaPI/PIext/nomhistadapter.cc

    r344 r361  
    298298}
    299299
    300 
     300//-------------------------------------------------------------------------
     301// Class Adaptateur d'objet (Pour NamedObjMgr) d'objet XNTuple
     302//-------------------------------------------------------------------------
     303
     304/* --Methode-- */
     305NOMAdapter_XNTuple::NOMAdapter_XNTuple(XNTuple* o)
     306  : NObjMgrAdapter(o)
     307{
     308mNt = o;
     309}
     310
     311/* --Methode-- */
     312NOMAdapter_XNTuple::~NOMAdapter_XNTuple()
     313{
     314}
     315
     316/* --Methode-- */
     317NObjMgrAdapter* NOMAdapter_XNTuple::Clone(AnyDataObj* o)
     318{
     319XNTuple* nt = dynamic_cast<XNTuple *>(o);
     320if (nt) return ( new NOMAdapter_XNTuple(nt) );
     321return ( new NObjMgrAdapter(o) );
     322}
     323
     324
     325/* --Methode-- */
     326void NOMAdapter_XNTuple::SavePPF(POutPersist& pos, string const & nom)
     327{
     328#ifdef SANS_EVOLPLANCK
     329// PEIDA-EROS L'histo est lui-meme PPersist
     330string tag = nom;  // A cause de const
     331mNt->Write(pos,0,tag);
     332#else
     333string s = typeid(*mObj).name();
     334cout << "NOMAdapter_XNTuple::SavePPF() - Error : Not supported for " << s << endl;
     335#endif
     336}
     337
     338/* --Methode-- */
     339void NOMAdapter_XNTuple::Print(ostream& os)
     340{
     341// os << mNt->Info();
     342mNt->Show(os);
     343}
     344
     345
     346/* --Methode-- */
     347NTupleInterface* NOMAdapter_XNTuple::GetNTupleInterface(bool& adel)
     348{
     349adel = false;
     350return(mNt);
     351}
     352
     353
  • trunk/SophyaPI/PIext/nomhistadapter.h

    r344 r361  
    1212#include "hisprof.h"
    1313#include "ntuple.h"
     14#include "xntuple.h"
    1415
    1516//-------------------------------------------------------------------------
     
    9596class NOMAdapter_NTuple : public NObjMgrAdapter {
    9697public:
    97                                 NOMAdapter_NTuple(NTuple* h = NULL);
     98                                NOMAdapter_NTuple(NTuple* nt = NULL);
    9899  virtual                       ~NOMAdapter_NTuple();
    99100
     
    111112};
    112113
     114//-------------------------------------------------------------------------
     115// Class Adaptateur d'objet (Pour NamedObjMgr) d'objet XNTuple
     116//-------------------------------------------------------------------------
     117
     118class NOMAdapter_XNTuple : public NObjMgrAdapter {
     119public:
     120                                NOMAdapter_XNTuple(XNTuple* nt = NULL);
     121  virtual                       ~NOMAdapter_XNTuple();
     122
     123  virtual NObjMgrAdapter*       Clone(AnyDataObj* o);
     124
     125  //  virtual void                      ReadFits(string const & flnm);
     126  //  virtual void                      SaveFits(string const & flnm);
     127  virtual void                  SavePPF(POutPersist& s, string const & nom);
     128
     129  virtual void                  Print(ostream& os);
     130  virtual NTupleInterface*      GetNTupleInterface(bool& adel);
     131
     132protected:
     133  XNTuple* mNt;
     134};
     135
    113136
    114137
  • trunk/SophyaPI/PIext/piacmd.cc

    r357 r361  
    1313#include "pistdimgapp.h"
    1414#include "nobjmgr.h"
     15#include "piafitting.h"
    1516
    1617#include PISTDWDG_H
     
    309310
    310311basexec = new PIABaseExecutor(this, omg, app);
     312fitexec = new PIAFitter(this, app);
    311313AddInterpreter(this);
    312314curcmdi = this;
     
    327329delete helpwin;
    328330if (curpiacmd == this)  curpiacmd = NULL;
     331delete basexec;
     332delete fitexec;
    329333}
    330334
  • trunk/SophyaPI/PIext/piacmd.h

    r351 r361  
    9797  CmdInterpreter* curcmdi;
    9898  CmdExecutor* basexec;
     99  CmdExecutor* fitexec;
    99100
    100101// Pour enregistrer la liste de commandes et leurs executeurs et le help
  • trunk/SophyaPI/PIext/piinit.cc

    r339 r361  
    2323  serv->RegisterClass(new Histo2D, new NOMAdapter_Histo2D );
    2424  serv->RegisterClass(new NTuple, new NOMAdapter_NTuple );
     25  serv->RegisterClass(new XNTuple, new NOMAdapter_XNTuple );
    2526
    2627  serv->RegisterClass(new GeneralFitData, new NOMAdapter_GeneralFitData );
  • trunk/SophyaPI/PIext/servnobjm.cc

    r357 r361  
    961961delete nt;
    962962mOmg->AddObj(gfd, nomgfd);
    963 return;
    964 }
    965 
    966 
    967 ///////////////////// Fit 1D et 2D //////////////////////////
    968 /* --Function static propre aux routines de fit 1D et 2D-- cmv 13/10/98 */
    969 struct DFOptions {
    970   int okres, okfun;
    971   int polcx,polcy; double xc,yc;
    972   double err_e, err_E;
    973   double stc2;
    974   int nstep;
    975   int lp,lpg;
    976   int i1,i2,j1,j2;
    977 };
    978 typedef struct DFOptions DFOPTIONS;
    979 static void DecodeFitsOptions(string par,string step,string min,string max,string opt
    980                       ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O);
    981 void DecodeFitsOptions(string par,string step,string min,string max,string opt
    982                       ,Vector& Par,Vector& Step,Vector& Min,Vector& Max,DFOPTIONS& O)
    983 //| Pour decoder les "string" et remplir les vecteurs du fit (cf commentaires dans Fit1D)
    984 {
    985 // set des vecteurs et decodage des string correspondantes
    986 int NParMax = 100;
    987 Par.Realloc(NParMax); Step.Realloc(NParMax);
    988 Min.Realloc(NParMax); Max.Realloc(NParMax);
    989 {
    990   Vector* v=NULL; string* s=NULL;
    991   {for(int i=0;i<NParMax;i++) {Par(i)=0.; Step(i)=1.; Min(i)=1.; Max(i)=-1.;}}
    992   for(int j=0;j<4;j++) {
    993     if(j==0)      {v=&Par; s=&par;}
    994     else if(j==1) {v=&Step; s=&step;}
    995     else if(j==2) {v=&Min; s=&min;}
    996     else if(j==3) {v=&Max; s=&max;}
    997     if(s->length()>0) *s += ",";
    998     for(int i=0;i<NParMax;i++) {
    999       if(s->length()<=0) break;
    1000       sscanf(s->c_str(),"%lf",&(*v)(i));
    1001       size_t p = s->find_first_of(',') + 1;
    1002       if(p>=s->length()) *s = ""; else *s = s->substr(p);
    1003     }
    1004   }
    1005 }
    1006 
    1007 // Decodage de options de opt
    1008 O.okres = O.okfun = 0;
    1009 O.polcx = O.polcy = 0;
    1010 O.xc = O.yc = 0.;
    1011 O.stc2 = 1.e-3;
    1012 O.nstep = 100;
    1013 O.err_e = O.err_E = -1.;
    1014 O.lp = 1; O.lpg = 0;
    1015 O.i1 = O.j1 = O.i2 = O.j2 = -1;
    1016 
    1017 if(opt.length()<=0) return;
    1018 opt = "," + opt + ",";
    1019 
    1020 if(strstr(opt.c_str(),",r,")) O.okres = 1;  // residus
    1021 if(strstr(opt.c_str(),",f,")) O.okfun = 1;  // fonction fittee
    1022 if(strstr(opt.c_str(),",x")) { // Demande de centrage (fit X=x-xc)
    1023   O.polcx = 2; // Le centrage est calcule automatiquement
    1024   size_t p = opt.find(",x");
    1025   size_t q = opt.find_first_of(',',p+1);
    1026   string dum = opt.substr(p,q-p);
    1027   if(dum.length()>2) {
    1028     sscanf(dum.c_str(),",x%lf",&O.xc);
    1029     O.polcx = 1; // Le centrage est fixe par la valeur lue
    1030   }
    1031 }
    1032 if(strstr(opt.c_str(),",y")) { // Demande de centrage (fit Y=y-yc)
    1033   O.polcy = 2; // Le centrage est calcule automatiquement
    1034   size_t p = opt.find(",y");
    1035   size_t q = opt.find_first_of(',',p+1);
    1036   string dum = opt.substr(p,q-p);
    1037   if(dum.length()>2) {
    1038     sscanf(dum.c_str(),",y%lf",&O.yc);
    1039     O.polcy = 1; // Le centrage est fixe par la valeur lue
    1040   }
    1041 }
    1042 if(strstr(opt.c_str(),",E")) { // Erreurs imposees a "sqrt(val)" ou "aa.b*sqrt(val)"
    1043   size_t p = opt.find(",E");
    1044   size_t q = opt.find_first_of(',',p+1);
    1045   string dum = opt.substr(p,q-p);
    1046   if(dum.length()>2) sscanf(dum.c_str(),",E%lf",&O.err_E);
    1047   if(O.err_E<=0.) O.err_E = 1.;
    1048   O.err_e=-1.;
    1049 }
    1050 if(strstr(opt.c_str(),",e")) { // Erreurs imposees a "1" ou "aa.b"
    1051   size_t p = opt.find(",e");
    1052   size_t q = opt.find_first_of(',',p+1);
    1053   string dum = opt.substr(p,q-p);
    1054   if(dum.length()>2) sscanf(dum.c_str(),",e%lf",&O.err_e);
    1055   if(O.err_e<=0.) O.err_e = 1.;
    1056   O.err_E=-1.;
    1057 }
    1058 if(strstr(opt.c_str(),",X")) { // Valeur du StopChi2
    1059   size_t p = opt.find(",X");
    1060   size_t q = opt.find_first_of(',',p+1);
    1061   string dum = opt.substr(p,q-p);
    1062   if(dum.length()>2) sscanf(dum.c_str(),",X%lf",&O.stc2);
    1063   if(O.stc2<=0.) O.stc2 = 1.e-3;
    1064 }
    1065 if(strstr(opt.c_str(),",N")) { // Nombre maximum d'iterations
    1066   size_t p = opt.find(",N");
    1067   size_t q = opt.find_first_of(',',p+1);
    1068   string dum = opt.substr(p,q-p);
    1069   if(dum.length()>2) sscanf(dum.c_str(),",N%d",&O.nstep);
    1070   if(O.nstep<2) O.nstep = 100;
    1071 }
    1072 if(strstr(opt.c_str(),",l")) { // niveau de print
    1073   size_t p = opt.find(",l");
    1074   size_t q = opt.find_first_of(',',p+1);
    1075   string dum = opt.substr(p,q-p);
    1076   float ab;
    1077   if(dum.length()>2) sscanf(dum.c_str(),",l%f",&ab);
    1078   if(ab<0) ab = 0.;
    1079   O.lp = (int) ab; O.lpg = int(10.*(ab-(float)O.lp+0.01));
    1080 }
    1081 if(strstr(opt.c_str(),",I")) { // intervalle de fit selon X
    1082   size_t p = opt.find(",I");
    1083   size_t q = opt.find_first_of(',',p+1);
    1084   string dum = opt.substr(p,q-p);
    1085   if(dum.length()>2) sscanf(dum.c_str(),",I%d/%d",&O.i1,&O.i2);
    1086 }
    1087 if(strstr(opt.c_str(),",J")) { // intervalle de fit selon Y
    1088   size_t p = opt.find(",J");
    1089   size_t q = opt.find_first_of(',',p+1);
    1090   string dum = opt.substr(p,q-p);
    1091   if(dum.length()>2) sscanf(dum.c_str(),",J%d/%d",&O.j1,&O.j2);
    1092 }
    1093 return;
    1094 }
    1095 
    1096 /* --Methode-- cmv 13/10/98 */
    1097 void  Services2NObjMgr::Fit12D(string& nom, string& func,
    1098                                string par,string step,string min,string max,
    1099                                string opt)
    1100 //| --------------- Fit d'objets a 1 et 2 dimensions ---------------
    1101 //| nom  : nom de l'objet qui peut etre:
    1102 //|        fit-1D:  Vector,Histo1D,HProf ou GeneraFitData(1D)
    1103 //|        fit-2D:  Matrix,Histo2D,Image<T> ou GeneraFitData(2D)
    1104 //| func : pnn : fit polynome degre nn avec classe Poly (lineaire) 1D ou 2D
    1105 //|      : Pnn : fit polynome degre nn avec GeneralFit (non-lineaire) 1D ou 2D
    1106 //|      : gnn : fit gaussienne (hauteur) + polynome de degre nn 1D
    1107 //|      : g   : fit gaussienne (hauteur) 1D
    1108 //|      : enn : fit exponentielle + polynome de degre nn 1D
    1109 //|      : e   : fit exponentielle 1D
    1110 //|      : Gnn : fit gaussienne (volume) + polynome de degre nn 1D
    1111 //|      : G   : fit gaussienne (volume) 1D
    1112 //|      :     : fit gaussienne+fond (volume) 2D
    1113 //|      : Gi  : fit gaussienne+fond integree (volume) 2D
    1114 //|      : d   : fit DL de gaussienne+fond (volume) 2D
    1115 //|      : di  : fit DL de gaussienne+fond integree (volume) 2D
    1116 //|      : D   : fit DL de gaussienne+fond avec coeff variable p6 (volume) 2D
    1117 //|      : Di  : fit DL de gaussienne+fond integree avec coeff variable p6 (volume) 2D
    1118 //|      : M   : fit Moffat+fond (expos=p6) (volume) 2D
    1119 //|      : Mi  : fit Moffat+fond integree (expos=p6) (volume) 2D
    1120 //| par  : p1,...,pn : valeur d'initialisation des parametres (def=0)
    1121 //| step : s1,...,sn : valeur des steps de depart (def=1)
    1122 //| min  : m1,...,mn : valeur des minima (def=1)
    1123 //| max  : M1,...,Mn : valeur des maxima (def=-1) (max<=min : pas de limite)
    1124 //| opt  : options "Eaa.b,eaa.b,f,r,caa.b,Xaa.b"
    1125 //|      f : generation d'un Objet identique contenant la fonction fittee
    1126 //|      r : generation d'un Objet identique contenant les residus
    1127 //|      Xaa.b : aa.b valeur du DXi2 d'arret (def=1.e-3)
    1128 //|      Naa : aa nombre maximum d'iterations (def=100)
    1129 //|      la.b : niveau "a.b" de print: a=niveau de print Fit1/2D
    1130 //|                                    b=niveau de debug GeneralFit
    1131 //|      Ii1/i2 numeros des bins X de l'histos utilises pour le fit [i1,i2]
    1132 //|2D    Jj1/j2 numeros des bins Y de l'histos utilises pour le fit [j1,j2]
    1133 //|      - L'erreur est celle associee a l'objet (si elle existe),
    1134 //|        elle est mise a 1 sinon, sauf si E... ou e... est precise:
    1135 //|      Eaa.b : si |val|>=1 erreur = aa.b*sqrt(|val|)
    1136 //|              si |val|<1  erreur = aa.b
    1137 //|              si aa.b <=0 alors aa.b=1.0
    1138 //|              E seul est equivalent a E1.0
    1139 //|      eaa.b : erreur = aa.b
    1140 //|              si aa.b <=0 alors aa.b=1.0
    1141 //|              e seul est equivalent a e1.0
    1142 //|      xaa.b : demande de centrage: on fit x-aa.b au lieu de x)
    1143 //|      x : demande de centrage: on fit x-xc au lieu de x
    1144 //|          avec xc=abscisse du milieu de l'histogramme
    1145 //|          Actif pour exp+poly 1D, poly 1D
    1146 //|                pour gauss+poly 1D, xc est le centre de la gaussienne.
    1147 //|2D    yaa.b et y : idem "xaa.b et x" mais pour y
    1148 {
    1149 AnyDataObj* obj=mOmg->GetObj(nom);
    1150 if (obj == NULL) {
    1151   cout<<"Services2NObjMgr::Fit12D() Error , Pas d'objet de nom "<<nom<<endl;
    1152   return;
    1153 }
    1154 if (!mImgapp)  return;
    1155 if(func.length()<=0)
    1156   {cout<<"Services2NObjMgr::Fit12D() Donnez un nom de fonction a fitter."<<endl;
    1157    return;}
    1158 string ctyp = typeid(*obj).name();
    1159 
    1160 int ndim = 0, nbinx=0, nbiny=0, ndata = 0;
    1161 Vector* v = NULL; Histo* h = NULL;
    1162 Matrix* m = NULL; Histo2D* h2 = NULL; RzImage* im = NULL;
    1163 GeneralFitData* g = NULL;
    1164 
    1165   // 1D
    1166 if (typeid(*obj) == typeid(Vector)) {
    1167   ndim = 1;
    1168   v = (Vector*) obj; nbinx = v->NElts(); nbiny = 1;
    1169   }
    1170 else if ( (typeid(*obj) == typeid(HProf)) || (typeid(*obj) == typeid(Histo)) ) {
    1171   ndim = 1;
    1172   h = (Histo*)  obj; nbinx = h->NBins(); nbiny = 1;
    1173   }
    1174 else if (typeid(*obj) == typeid(Matrix)) {
    1175   ndim = 2;
    1176   m = (Matrix*) obj; nbinx = m->NCol(); nbiny = m->NRows();
    1177   }
    1178 else if (typeid(*obj) == typeid(Histo2D)) {
    1179   ndim = 2;
    1180   h2 = (Histo2D*) obj; nbinx = h2->NBinX(); nbiny = h2->NBinY();
    1181   }
    1182 else if (typeid(*obj) == typeid(GeneralFitData)) {
    1183   g = (GeneralFitData*) obj; nbinx = g->NData(); nbiny = 1;
    1184   if(     g->NVar()==1) ndim = 1;
    1185   else if(g->NVar()==2) ndim = 2;
    1186   else {
    1187      cout<<"GeneralFitData ne peut avoir que 1 ou 2 variables d'abscisse: "
    1188          <<((GeneralFitData*) obj)->NVar()<<endl; return; }
    1189   }
    1190 else if (dynamic_cast<RzImage*>(obj)) {
    1191   ndim = 2;
    1192   im = (RzImage*) obj; nbinx = im->XSize(); nbiny = im->YSize();
    1193   }
    1194 else  {
    1195   cout<<"Services2NObjMgr::Fit12D() Error , Objet n'est pas un "
    1196       <<"Histo1D/HProf/Vector/Histo2D/Image/Matrix/GeneralFitData "<<ctyp<<endl;
    1197   return;
    1198   }
    1199 
    1200 ndata = nbinx*nbiny;
    1201 if(ndata<=0)
    1202   {cout<<"L'objet a "<<nbinx<<","<<nbiny<<" bins ("<<ndata<<")"<<endl; return;}
    1203 
    1204 // Decodage des options et des parametres, mise en forme
    1205 Vector Par(1); Vector Step(1); Vector Min(1); Vector Max(1); DFOPTIONS O;
    1206 DecodeFitsOptions(par,step,min,max,opt,Par,Step,Min,Max,O);
    1207 O.i1 = (O.i1<0||O.i1>=nbinx)? 0: O.i1;
    1208 O.i2 = (O.i2<0||O.i2>=nbinx||O.i2<O.i1)? nbinx-1: O.i2;
    1209 if(ndim>=2) {
    1210   O.j1 = (O.j1<0||O.j1>=nbiny)? 0: O.j1;
    1211   O.j2 = (O.j2<0||O.j2>=nbiny||O.j2<O.j1)? nbiny-1: O.j2;
    1212 } else O.j2 = O.j1 = 0;
    1213 if(O.polcx==2) {
    1214   if(v||m)    O.xc = (O.i2-O.i1+1)/2.;
    1215   else if(h)  O.xc = (h->XMin()+h->XMax())/2.;
    1216   else if(h2) O.xc = (h2->XMin()+h2->XMax())/2.;
    1217   else if(g)  {double mini,maxi; g->GetMinMax(2,mini,maxi); O.xc=(mini+maxi)/2.;}
    1218   else if(im) {O.xc = im->XOrg() * im->XPxSize()*(O.i2-O.i1+1)/2.;}
    1219 }
    1220 if(O.polcy==2 && ndim>=2) {
    1221   if(m)  O.yc = (O.j2-O.j1+1)/2.;
    1222   if(h2) O.yc = (h2->YMin()+h2->YMax())/2.;
    1223   if(g)  {double mini,maxi; g->GetMinMax(12,mini,maxi); O.yc=(mini+maxi)/2.;}
    1224   if(im) {O.yc = im->YOrg() * im->YPxSize()*(O.j2-O.j1+1)/2.;}
    1225 }
    1226 if(O.lp>0)
    1227   cout<<"Fit["<<nbinx<<","<<nbiny<<"] ("<<ndata<<") dim="<<ndim<<":"
    1228       <<" Int=["<<O.i1<<","<<O.i2<<"],["<<O.j1<<","<<O.j2<<"]"<<endl
    1229       <<" Cent="<<O.polcx<<","<<O.polcy<<","<<O.xc<<"+x"<<","<<O.yc<<"+y"
    1230       <<" TypE="<<O.err_e<<","<<O.err_E
    1231       <<" StpX2="<<O.stc2<<" Nstep="<<O.nstep
    1232       <<" lp,lpg="<<O.lp<<","<<O.lpg<<endl;
    1233 
    1234 ///////////////////////////////////
    1235 // Remplissage de GeneralFitData //
    1236 ///////////////////////////////////
    1237 GeneralFitData mydata(ndim,ndata,0);
    1238 {for(int i=O.i1;i<=O.i2;i++) for(int j=O.j1;j<=O.j2;j++) {
    1239   double x,y,f,e;
    1240 
    1241   if(v)
    1242     {x= (double) i; f=(*v)(i); e=1.;}
    1243   else if(h)
    1244     {x=h->BinCenter(i); f=(*h)(i); e=(h->HasErrors())?h->Error(i):1.;}
    1245   else if(m)
    1246     {x=(double) i; y=(double) j; f=(*m)(j,i); e=1.;}
    1247   else if(h2)
    1248     {float xf,yf; h2->BinCenter(i,j,xf,yf); x=(double)xf; y=(double)yf;
    1249      f=(*h2)(i,j); e=(h2->HasErrors())?h2->Error(i,j):1.;}
    1250   else if(im)
    1251     {x=im->XOrg()+(i+0.5)*im->XPxSize(); y=im->YOrg()+(j+0.5)*im->YPxSize();
    1252      f=im->DValue(i,j); e=1.;}
    1253   else if(g&&ndim==1) {x= g->X(i); f=g->Val(i); e=g->EVal(i);}
    1254   else if(g&&ndim==2) {x= g->X(i); y= g->Y(i); f=g->Val(i); e=g->EVal(i);}
    1255   else x=y=f=e=0.;
    1256 
    1257   // Gestion des erreurs a utiliser
    1258   if(O.err_e>0.) e=O.err_e;
    1259   else if(O.err_E>0.) {e=(f<-1.||f>1.)?O.err_E*sqrt(fabs(f)):O.err_E;}
    1260 
    1261   // Remplissage de generalfit
    1262   if(func[0]=='p') {x -= O.xc; if(ndim>=2) y -= O.yc;}
    1263   if(ndim==1)      mydata.AddData1(x,f,e);
    1264   else if(ndim==2) mydata.AddData2(x,y,f,e);
    1265 }}
    1266 if(mydata.NData()<=0)
    1267   {cout<<"Pas de donnees dans GeneralFitData: "<<mydata.NData()<<endl;
    1268    return;}
    1269 if(O.lpg>1) {
    1270   mydata.PrintStatus();
    1271   mydata.PrintData(0);
    1272   mydata.PrintData(mydata.NData()-1);
    1273 }
    1274 
    1275 ////////////////////////////////////////////
    1276 // Identification de la fonction a fitter //
    1277 ////////////////////////////////////////////
    1278 GeneralFunction* myfunc = NULL;
    1279 if(func[0]=='p' && ndim==1) {
    1280   // Fit de polynome sans passer par les GeneralFit
    1281   int degre = 0;
    1282   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1283   cout<<"Fit (lineaire) 1D polynome de degre "<<degre<<endl;
    1284   Poly p1(0);
    1285   double c2rl = mydata.PolFit(0,p1,degre);
    1286   cout<<"C2r_lineaire = "<<c2rl<<endl;
    1287   if(O.lp>0) cout<<p1<<endl;
    1288   return;
    1289 
    1290 } else if(func[0]=='P' && ndim==1) {
    1291   // Fit de polynome
    1292   int degre = 0;
    1293   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1294   cout<<"Fit polynome 1D de degre "<<degre<<endl;
    1295   Polyn1D* myf = new Polyn1D(degre,O.xc);
    1296   myfunc = myf;
    1297 
    1298 } else if(func[0]=='e' && ndim==1) {
    1299   // Fit d'exponentielle
    1300   int degre =-1;
    1301   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1302   cout<<"Fit d'exponentielle+polynome 1D de degre "<<degre<<endl;
    1303   Exp1DPol* myf;
    1304   if(degre>=0) myf = new Exp1DPol((unsigned int)degre,O.xc);
    1305        else    myf = new Exp1DPol(O.xc);
    1306   myfunc = myf;
    1307 
    1308 } else if(func[0]=='g' && ndim==1) {
    1309   // Fit de gaussienne en hauteur
    1310   int degre =-1;
    1311   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1312   cout<<"Fit de Gaussienne_en_hauteur+polynome 1D de degre "<<degre<<endl;
    1313   Gauss1DPol* myf;
    1314   if(degre>=0) myf = new Gauss1DPol((unsigned int)degre,((O.polcx)?true:false));
    1315   else { bool bfg = (O.polcx)?true:false;   myf = new Gauss1DPol(bfg); }
    1316   myfunc = myf;
    1317 
    1318 } else if(func[0]=='G' && ndim==1) {
    1319   // Fit de gaussienne en volume
    1320   int degre =-1;
    1321   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1322   cout<<"Fit de Gaussienne_en_volume+polynome 1D de degre "<<degre<<endl;
    1323   GaussN1DPol* myf;
    1324   if(degre>=0) myf = new GaussN1DPol((unsigned int)degre,((O.polcx)?true:false));
    1325   else  { bool bfg = (O.polcx)?true:false;   myf = new GaussN1DPol(bfg); }
    1326   myfunc = myf;
    1327 
    1328 } else if(func[0]=='p' && ndim==2) {
    1329   // Fit de polynome 2D sans passer par les GeneralFit
    1330   int degre = 0;
    1331   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1332   cout<<"Fit (lineaire) polynome 2D de degre "<<degre<<endl;
    1333   Poly2 p2(0);
    1334   double c2rl = mydata.PolFit(0,1,p2,degre);
    1335   cout<<"C2r_lineaire = "<<c2rl<<endl;
    1336   if(O.lp>0) cout<<p2<<endl;
    1337   return;
    1338 
    1339 } else if(func[0]=='P' && ndim==2) {
    1340   // Fit de polynome 2D
    1341   int degre = 0;
    1342   if(func.length()>1) sscanf(func.c_str()+1,"%d",&degre);
    1343   cout<<"Fit polynome 2D de degre "<<degre<<endl;
    1344   Polyn2D* myf = new Polyn2D(degre,O.xc,O.yc);
    1345   myfunc = myf;
    1346 
    1347 } else if(func[0]=='G' && ndim==2) {
    1348   // Fit de gaussienne+fond en volume
    1349   int integ = 0;
    1350   if(func.length()>1) if(func[1]=='i') integ=1;
    1351   cout<<"Fit de Gaussienne+Fond 2D integ="<<integ<<endl;
    1352   if(integ) {GauRhInt2D* myf = new GauRhInt2D; myfunc = myf;}
    1353     else    {GauRho2D* myf = new GauRho2D; myfunc = myf;}
    1354 
    1355 } else if(func[0]=='d' && ndim==2) {
    1356   // Fit de DL gaussienne+fond en volume
    1357   int integ = 0;
    1358   if(func.length()>1) if(func[1]=='i') integ=1;
    1359   cout<<"Fit de DL de Gaussienne+Fond 2D integ="<<integ<<endl;
    1360   if(integ) {GdlRhInt2D* myf = new GdlRhInt2D; myfunc = myf;}
    1361     else    {GdlRho2D* myf = new GdlRho2D; myfunc = myf;}
    1362 
    1363 } else if(func[0]=='D' && ndim==2) {
    1364   // Fit de DL gaussienne+fond avec coeff variable p6 en volume
    1365   int integ = 0;
    1366   if(func.length()>1) if(func[1]=='i') integ=1;
    1367   cout<<"Fit de DL de Gaussienne+Fond avec coeff variable (p6) 2D integ="<<integ<<endl;
    1368   if(integ) {Gdl1RhInt2D* myf = new Gdl1RhInt2D; myfunc = myf;}
    1369     else    {Gdl1Rho2D* myf = new Gdl1Rho2D; myfunc = myf;}
    1370 
    1371 } else if(func[0]=='M' && ndim==2) {
    1372   // Fit de Moffat+fond (volume)
    1373   int integ = 0;
    1374   if(func.length()>1) if(func[1]=='i') integ=1;
    1375   cout<<"Fit de Moffat+Fond (expos=p6) 2D integ="<<integ<<endl;
    1376   if(integ) {MofRhInt2D* myf = new MofRhInt2D; myfunc = myf;}
    1377     else    {MofRho2D* myf = new MofRho2D; myfunc = myf;}
    1378 
    1379 } else {
    1380   cout<<"Fonction "<<func<<" inconnue pour la dim "<<ndim<<endl;
    1381   return;
    1382 }
    1383 
    1384 /////////////////////////
    1385 // Fit avec generalfit //
    1386 /////////////////////////
    1387 if(myfunc->NPar()>Par.NElts())
    1388   {cout<<"Trop de parametres: "<<myfunc->NPar()<<">"<<Par.NElts()<<endl;
    1389   if(myfunc) delete myfunc; return;}
    1390 GeneralFit myfit(myfunc);
    1391 myfit.SetDebug(O.lpg);
    1392 myfit.SetData(&mydata);
    1393 myfit.SetStopChi2(O.stc2);
    1394 myfit.SetMaxStep(O.nstep);
    1395 {for(int i=0;i<myfunc->NPar();i++) {
    1396   char str[10];
    1397   sprintf(str,"P%d",i);
    1398   myfit.SetParam(i,str,Par(i),Step(i),Min(i),Max(i));
    1399 }}
    1400 if(O.lp>1) myfit.PrintFit();
    1401 double c2r = -1.;
    1402 int rcfit = (double) myfit.Fit();
    1403 if(O.lp>0) myfit.PrintFit();
    1404 if(rcfit>0) {
    1405   c2r = myfit.GetChi2Red();
    1406   cout<<"C2r_Reduit = "<<c2r<<" nstep="<<myfit.GetNStep()<<" rc="<<rcfit<<endl;
    1407   Vector ParFit(myfunc->NPar());
    1408   for(int i=0;i<myfunc->NPar();i++) ParFit(i)=myfit.GetParm(i);
    1409 } else {
    1410   cout<<"echec Fit, rc = "<<rcfit<<"  nstep="<<myfit.GetNStep()<<endl;
    1411   myfit.PrintFitErr(rcfit);
    1412 }
    1413 
    1414 // Mise a disposition des resultats
    1415 if(rcfit>=0 && myfunc && (O.okres>0||O.okfun>0)) {
    1416   string nomres = nom + "res";
    1417   string nomfun = nom + "fun";
    1418   if(v) {
    1419     if(O.okres) {Vector* ob = v->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1420     if(O.okfun) {Vector* ob = v->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1421   } else if(h) {
    1422     if(O.okres) {Histo* ob = h->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1423     if(O.okfun) {Histo* ob = h->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1424   } else if(m) {
    1425     if(O.okres) {Matrix* ob = m->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1426     if(O.okfun) {Matrix* ob = m->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1427   } else if(h2) {
    1428     if(O.okres) {Histo2D* ob = h2->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1429     if(O.okfun) {Histo2D* ob = h2->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1430   } else if(im) {
    1431     if(O.okres) {RzImage* ob = im->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1432     if(O.okfun) {RzImage* ob = im->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1433   } else if(g) {
    1434     if(O.okres) {GeneralFitData* ob = g->FitResidus(myfit);  if(ob) mOmg->AddObj(ob,nomres);}
    1435     if(O.okfun) {GeneralFitData* ob = g->FitFunction(myfit); if(ob) mOmg->AddObj(ob,nomfun);}
    1436   }
    1437 }
    1438 
    1439 // Nettoyage
    1440 if(myfunc) delete myfunc;
    1441963return;
    1442964}
     
    19321454  case ClassId_NTuple :
    19331455    return("NTuple");
     1456  case ClassId_XNTuple :
     1457    return("XNTuple");
    19341458  case ClassId_GeneralFitData :
    19351459    return("GeneralFitData");
  • trunk/SophyaPI/PIext/servnobjm.h

    r357 r361  
    9191                                                 string const & funcname);
    9292
    93 //   Methodes de fit -  CMV , deux methodes H1 H2 ou 1 seule ??) ...
    94   virtual void          Fit12D(string & nom, string& func,
    95                               string par,string step,string min,string max,string opt);
    96 
    9793// Calcul d'expressions d'interface NTuple pour les objets
    9894  void          ComputeExpressions(NObjMgrAdapter* obja, string& expx, string& expy,
Note: See TracChangeset for help on using the changeset viewer.