Index: /trunk/Cosmo/SimLSS/cmvgsm2cube.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvgsm2cube.cc	(revision 3982)
+++ /trunk/Cosmo/SimLSS/cmvgsm2cube.cc	(revision 3983)
@@ -11,6 +11,7 @@
 #include "skymap.h"
 
-// set fname = ( `cat gsm820948.list | awk '{print $3}' | sed 's/data/ppf/'` )
-// cmvgsm2cube -o gsm820948_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
+// set rname = gsm820948
+// set fname = ( `cat ${rname}.list | awk '{print $3}' | sed 's/data/ppf/'` )
+// cmvgsm2cube -o ${rname}_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
 
 void usage(void)
@@ -19,6 +20,6 @@
 <<"cmvgsm2cube [options] liste_des_spheres_healpix"<<endl
 <<" -o outname.ppf : ppf cube output file name (def=cmvgsm2cube.ppf)"<<endl
-<<" -p nphi,dphi,phi0 : phi0 -> phi0+nphi*dphi (degres)"<<endl
-<<" -t ntheta,dtheta,theta0 : theta0 -> theta0+ntheta*dtheta (degres)"<<endl
+<<" -p nphi,dphi,phi0 : phi0 -> phi0 + nphi*dphi (degres)"<<endl
+<<" -t ntheta,dtheta,theta0 : theta0 -> theta0 + ntheta*dtheta (degres)"<<endl
 <<endl;
 }
@@ -66,29 +67,57 @@
  cout<<"Number of spheres "<<sphname.size()<<endl;
 
- TArray<r_4> Cube(sphname.size(),nphi,ntheta);
- Cube = 0.;
+ TVector<r_4> Theta(ntheta);
+ int nthgood = 0;
+ float thmin = 1.e30, thmax = -1.e30;
+ for(int it=0;it<ntheta;it++) {
+   Theta(it) = -9999.;
+   double t = theta0 + it*dtheta; t *= M_PI/180.;
+   if(t<0. || t>M_PI) continue;  // [0,Pi]
+   Theta(it) = t;
+   nthgood += 1;
+   if(t<thmin) thmin = t;
+   if(t>thmax) thmax = t;
+ }
+ cout<<"Theta: ngood="<<nthgood<<" ["<<thmin*180/M_PI<<","<<thmax*180/M_PI
+     <<"] middle="<<0.5*(thmin+thmax)*180/M_PI<<endl;
+ 
+ TVector<r_4> Phi(nphi);
+ float phmin = 1.e30, phmax = -1.e30;
+ for(int ip=0;ip<nphi;ip++) {
+   double p = phi0 + ip*dphi; p *= M_PI/180.;
+   while(p<0.) p += 2.*M_PI;
+   while(p>=2.*M_PI) p -= 2.*M_PI;
+   Phi(ip) = p;
+   if(p<phmin) phmin = p;
+   if(p>phmax) phmax = p;
+ }
+ cout<<"Phi: ["<<phmin*180/M_PI<<","<<phmax*180/M_PI
+     <<"] middle="<<0.5*(phmin+phmax)*180/M_PI<<endl;
+
+
+ TVector<r_4> Freq(sphname.size()); Freq = -1.;
+ TArray<r_4> Cube(sphname.size(),nphi,ntheta); Cube = 0.;
  for(unsigned int is=0;is<sphname.size();is++) {
-
    cout<<"... "<<sphname[is]<<endl;
    PInPersist pis(sphname[is]);
    SphereHEALPix<r_4> sph;
    pis>>PPFNameTag("sph")>>sph;
-
-   for(int it=0;it<ntheta;it++) {
-     double t = theta0 + it*dtheta; t *= M_PI/180.;
-     if(t<0. || t>M_PI) continue;
-   for(int ip=0;ip<nphi;ip++) {
-     double p = phi0 + ip*dphi; p *= M_PI/180.;
-     while(p<0.) p += 2.*M_PI;
-     while(p>=2.*M_PI) p -= 2.*M_PI;
-     int_4 ksph = sph.PixIndexSph(t,p);
-     Cube(is,ip,it) = sph(ksph);
+   Freq(is) = (r_4)sph.Info()["FMHz"];
+   for(int it=0;it<Theta.Size();it++) {
+     double t = Theta(it);
+     if(t<0.) continue;
+     for(int ip=0;ip<Phi.Size();ip++) {
+       double p = Phi(ip);
+       int_4 ksph = sph.PixIndexSph(t,p);
+       Cube(is,ip,it) = sph(ksph);
+     }
    }
-   }
-
  }
  
  POutPersist pos(ppfrn);
  pos<<PPFNameTag("cube")<<Cube;
+ pos<<PPFNameTag("theta")<<Theta;
+ pos<<PPFNameTag("phi")<<Phi;
+ pos<<PPFNameTag("freq")<<Freq;
 
  return 0;
@@ -100,4 +129,8 @@
 print cube
 
+n/plot phi.val*180./M_PI%n ! ! "cpts"
+n/plot theta.val*180./M_PI%n ! ! "cpts"
+disp freq
+
 objaoper cube sliceyz 0
 objaoper cube sliceyz 128
Index: /trunk/Cosmo/SimLSS/cmvgsmcfit.cc
===================================================================
--- /trunk/Cosmo/SimLSS/cmvgsmcfit.cc	(revision 3982)
+++ /trunk/Cosmo/SimLSS/cmvgsmcfit.cc	(revision 3983)
@@ -24,14 +24,16 @@
 class MyFCN : public FCNBase {
 public:
-  MyFCN(double nu0,double dnu,TVector<r_8>& V,TVector<r_8>& EV);
+  MyFCN(TVector<r_4>& FreqMHz,TVector<r_4>& V,TVector<r_4>& EV);
   void InitFit(MnUserParameters& upar);
-  double Func(int ix,const vector<double>& par) const;
+  double Func(double fr,const vector<double>& par) const;
   double operator()(const vector<double>& par) const;
   double Up(void) const {return 1.;}
+  void SetNu0(double nu0Mhz) {_nu0 = nu0Mhz;}
   //protected:
+  double _nu0;
   mutable int _ndata;
-  double _nu0, _dnu;
-  TVector<r_8> _V;
-  TVector<r_8> _EV;
+  TVector<r_4> _F;
+  TVector<r_4> _V;
+  TVector<r_4> _EV;
 };
 } }  // namespace ROOT + Minuit2
@@ -40,9 +42,15 @@
 int main(int narg,char *arg[])
 {
-  double nu0 = 820.0, dnu = 0.5; // MHz
+  int nslicemax = -1;
   if(narg<=1) {
-    cout<<"cmvgsmcfit cubefile.ppf"<<endl;
+    cout<<"cmvgsmcfit cubefile.ppf [nslicemax]"<<endl;
     return -1;
   }
+
+  if(narg>2) nslicemax = atoi(arg[2]);
+  cout<<"Fitting "<<arg[1]<<" , nslicemax="<<nslicemax<<endl;
+
+int rc = 0;
+try {
 
   POutPersist pos("cmvgsmcfit.ppf");
@@ -50,8 +58,9 @@
   TArray<r_4> Cube;
   pis>>PPFNameTag("cube")>>Cube;
-  TVector<r_8> V(Cube.SizeX()), EV(Cube.SizeX());
-
-  const int nvar = 6;
-  const char *vname[nvar] = {"x2r","ok","ph","th","A","N"};
+  TVector<r_4> Freq;
+  pis>>PPFNameTag("freq")>>Freq;
+
+  const int nvar = 8;
+  const char *vname[nvar] = {"x2r","ok","ph","th","v0","vn","A","N"};
   NTuple nt(nvar,vname);
   double xnt[nvar]; for(int i=0;i<nvar;i++) xnt[i] = 0.;
@@ -59,6 +68,8 @@
   int nfit = 0;
   for(int j=0;j<Cube.SizeY();j++) {
-    cout<<"..."<<j<<endl;
+    cout<<"...j = "<<j<<endl;
   for(int k=0;k<Cube.SizeZ();k++) {
+
+    TVector<r_4> V(Cube.SizeX()), EV(Cube.SizeX());
     for(int i=0;i<Cube.SizeX();i++) {
       V(i) = Cube(i,j,k);
@@ -66,5 +77,6 @@
     }
 
-    MyFCN myfcn(nu0,dnu,V,EV);
+    MyFCN myfcn(Freq,V,EV);
+    myfcn.SetNu0(Freq(0));
     MnUserParameters upar;
     myfcn.InitFit(upar);
@@ -89,13 +101,14 @@
     xnt[0] = xi2r; xnt[1] = (okfit)? 1.: 0.;
     xnt[2] = j; xnt[3] = k;
-    for(unsigned int i=0;i<par.size();i++) xnt[4+i] = par[i];
+    xnt[4] = V(0); xnt[5] = V(V.Size()-1);
+    for(unsigned int i=0;i<par.size();i++) xnt[6+i] = par[i];
     nt.Fill(xnt);
     //vector<double> err = ustate.Errors();
-    for(int i=0;i<Cube.SizeX();i++) Cube(i,j,k) -= myfcn.Func(i,par);
+    for(int i=0;i<Cube.SizeX();i++) Cube(i,j,k) -= myfcn.Func(Freq(i),par);
 
     if(nfit<10) {
       TVector<r_8> F(Cube.SizeX()), R(Cube.SizeX());
       for(int i=0;i<Cube.SizeX();i++) {
-        F(i) = myfcn.Func(i,par);
+        F(i) = myfcn.Func(Freq(i),par);
         R(i) = V(i) - F(i);
       }
@@ -108,28 +121,46 @@
       pos<<PPFNameTag(str)<<R;
     }
+
     nfit++;
-    if(nfit>10) break;
+
   }
+    if(nslicemax>0 && j>=nslicemax) break;
   }
 
+  cout<<"nfit = "<<nfit<<endl;
   pos<<PPFNameTag("nt")<<nt;
   pos<<PPFNameTag("cublis")<<Cube;
 
+//-----------------------------------------------------------------
+} catch (PException& exc) {
+  cerr<<"cmvgsmcfit.cc catched PException"<<exc.Msg()<<endl;
+  rc = 77;
+} catch (std::exception& sex) {
+  cerr << "cmvgsmcfit.cc std::exception :" 
+       << (string)typeid(sex).name() << "\n msg= " << sex.what() << endl;
+  rc = 78;
+} catch (...) {
+  cerr << "cmvgsmcfit.cc catched unknown (...) exception  " << endl; 
+  rc = 79;
+}
+
+return rc;
 }
 
 //--------------------------------------------------
-MyFCN::MyFCN(double nu0, double dnu,TVector<r_8>& V,TVector<r_8>& EV)
-  : _V(V,false), _EV(EV,false)
-{
+MyFCN::MyFCN(TVector<r_4>& FreqMHz,TVector<r_4>& V,TVector<r_4>& EV)
+  : _F(FreqMHz,false), _V(V,false), _EV(EV,false)
+{
+  if(_F.Size()==0 || _F.Size()!=_V.Size() || _F.Size()!=_EV.Size())
+    throw SzMismatchError("MyFCN::MyFCN_Error: incompatible F / V / EV sizes");
+  _nu0 = _F(0);
   _ndata =_V.Size();
-  _nu0 = nu0;
-  _dnu = dnu;
 }
 
 void MyFCN::InitFit(MnUserParameters& upar)
 {
-  int n = _V.Size()-1;
+  int n = _F.Size() - 1;
   double a = _V(0);
-  double ni = log(_V(n)/a) / log(1.+n*_dnu/_nu0);
+  double ni = log(_V(n)/a) / log(_F(n)/_nu0);
 
   upar.Add("A",a,a/50.);
@@ -137,9 +168,8 @@
 }
 
-double MyFCN::Func(int ix,const vector<double>& par) const
-// f = A*(nu/nu0)^N   avec  nu = nu0 + i*dnu
-{
- double x = (_nu0 + ix*_dnu)/_nu0;
- double f = par[0]*pow(x,par[1]);
+double MyFCN::Func(double fr,const vector<double>& par) const
+// v = A*(nu/nu0)^N
+{
+ double f = par[0]*pow(fr/_nu0,par[1]);
  return f;
 }
@@ -153,5 +183,5 @@
    if(e2<=0.) continue;
    double v = _V(i);
-   double f = Func(i,par);
+   double f = Func(_F(i),par);
    xi2 += (v-f)*(v-f)/e2;
    _ndata++;
@@ -166,13 +196,26 @@
 set j 0
 set k 0
+zone 1 2
 n/plot v_${j}_${k}.val%n ! ! "nsta cpts"
 n/plot f_${j}_${k}.val%n ! ! "nsta cpts same red"
 n/plot r_${j}_${k}.val%n ! ! "nsta cpts"
 
-n/plot nt.ok%_nl ! ! "cpts"
-n/plot nt.x2r%_nl ! ! "cpts"
-
-n/plot nt.A%_nl ! ! "cpts"
-n/plot nt.N%_nl ! ! "cpts"
+zone 1 2
+n/plot nt.ok%_nl ! ! "crossmarker3"
+n/plot nt.x2r%_nl ! ! "crossmarker3"
+
+zone 1 2
+n/plot nt.v0%_nl ! ! "crossmarker3"
+n/plot nt.vn%_nl ! ! "plusmarker3 same red"
+n/plot nt.log10(v0/vn)%v0 ! ! "crossmarker3"
+
+
+zone
+n/plot nt.A%_nl ! ! "crossmarker3"
+n/plot nt.N%_nl ! ! "crossmarker3"
+
+n/plot nt.A%v0 ! ! "crossmarker3"
+n/plot nt.(A-v0)/v0%_nl ! ! "crossmarker3"
+n/plot nt.N%A ! ! "crossmarker3"
 
 set n 0
