Changeset 3983 in Sophya for trunk/Cosmo/SimLSS


Ignore:
Timestamp:
May 4, 2011, 7:13:16 PM (14 years ago)
Author:
cmv
Message:

amelioration des reconstructions et fits du GSM, cmv 4/5/2011

Location:
trunk/Cosmo/SimLSS
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/SimLSS/cmvgsm2cube.cc

    r3790 r3983  
    1111#include "skymap.h"
    1212
    13 // set fname = ( `cat gsm820948.list | awk '{print $3}' | sed 's/data/ppf/'` )
    14 // cmvgsm2cube -o gsm820948_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
     13// set rname = gsm820948
     14// set fname = ( `cat ${rname}.list | awk '{print $3}' | sed 's/data/ppf/'` )
     15// cmvgsm2cube -o ${rname}_cube.ppf -p 360,0.1666667,160. -t 360,0.1666667,30. $fname
    1516
    1617void usage(void)
    … …  
    1920<<"cmvgsm2cube [options] liste_des_spheres_healpix"<<endl
    2021<<" -o outname.ppf : ppf cube output file name (def=cmvgsm2cube.ppf)"<<endl
    21 <<" -p nphi,dphi,phi0 : phi0 -> phi0+nphi*dphi (degres)"<<endl
    22 <<" -t ntheta,dtheta,theta0 : theta0 -> theta0+ntheta*dtheta (degres)"<<endl
     22<<" -p nphi,dphi,phi0 : phi0 -> phi0 + nphi*dphi (degres)"<<endl
     23<<" -t ntheta,dtheta,theta0 : theta0 -> theta0 + ntheta*dtheta (degres)"<<endl
    2324<<endl;
    2425}
    … …  
    6667 cout<<"Number of spheres "<<sphname.size()<<endl;
    6768
    68  TArray<r_4> Cube(sphname.size(),nphi,ntheta);
    69  Cube = 0.;
     69 TVector<r_4> Theta(ntheta);
     70 int nthgood = 0;
     71 float thmin = 1.e30, thmax = -1.e30;
     72 for(int it=0;it<ntheta;it++) {
     73   Theta(it) = -9999.;
     74   double t = theta0 + it*dtheta; t *= M_PI/180.;
     75   if(t<0. || t>M_PI) continue;  // [0,Pi]
     76   Theta(it) = t;
     77   nthgood += 1;
     78   if(t<thmin) thmin = t;
     79   if(t>thmax) thmax = t;
     80 }
     81 cout<<"Theta: ngood="<<nthgood<<" ["<<thmin*180/M_PI<<","<<thmax*180/M_PI
     82     <<"] middle="<<0.5*(thmin+thmax)*180/M_PI<<endl;
     83 
     84 TVector<r_4> Phi(nphi);
     85 float phmin = 1.e30, phmax = -1.e30;
     86 for(int ip=0;ip<nphi;ip++) {
     87   double p = phi0 + ip*dphi; p *= M_PI/180.;
     88   while(p<0.) p += 2.*M_PI;
     89   while(p>=2.*M_PI) p -= 2.*M_PI;
     90   Phi(ip) = p;
     91   if(p<phmin) phmin = p;
     92   if(p>phmax) phmax = p;
     93 }
     94 cout<<"Phi: ["<<phmin*180/M_PI<<","<<phmax*180/M_PI
     95     <<"] middle="<<0.5*(phmin+phmax)*180/M_PI<<endl;
     96
     97
     98 TVector<r_4> Freq(sphname.size()); Freq = -1.;
     99 TArray<r_4> Cube(sphname.size(),nphi,ntheta); Cube = 0.;
    70100 for(unsigned int is=0;is<sphname.size();is++) {
    71 
    72101   cout<<"... "<<sphname[is]<<endl;
    73102   PInPersist pis(sphname[is]);
    74103   SphereHEALPix<r_4> sph;
    75104   pis>>PPFNameTag("sph")>>sph;
    76 
    77    for(int it=0;it<ntheta;it++) {
    78      double t = theta0 + it*dtheta; t *= M_PI/180.;
    79      if(t<0. || t>M_PI) continue;
    80    for(int ip=0;ip<nphi;ip++) {
    81      double p = phi0 + ip*dphi; p *= M_PI/180.;
    82      while(p<0.) p += 2.*M_PI;
    83      while(p>=2.*M_PI) p -= 2.*M_PI;
    84      int_4 ksph = sph.PixIndexSph(t,p);
    85      Cube(is,ip,it) = sph(ksph);
     105   Freq(is) = (r_4)sph.Info()["FMHz"];
     106   for(int it=0;it<Theta.Size();it++) {
     107     double t = Theta(it);
     108     if(t<0.) continue;
     109     for(int ip=0;ip<Phi.Size();ip++) {
     110       double p = Phi(ip);
     111       int_4 ksph = sph.PixIndexSph(t,p);
     112       Cube(is,ip,it) = sph(ksph);
     113     }
    86114   }
    87    }
    88 
    89115 }
    90116 
    91117 POutPersist pos(ppfrn);
    92118 pos<<PPFNameTag("cube")<<Cube;
     119 pos<<PPFNameTag("theta")<<Theta;
     120 pos<<PPFNameTag("phi")<<Phi;
     121 pos<<PPFNameTag("freq")<<Freq;
    93122
    94123 return 0;
    … …  
    100129print cube
    101130
     131n/plot phi.val*180./M_PI%n ! ! "cpts"
     132n/plot theta.val*180./M_PI%n ! ! "cpts"
     133disp freq
     134
    102135objaoper cube sliceyz 0
    103136objaoper cube sliceyz 128

  • trunk/Cosmo/SimLSS/cmvgsmcfit.cc

    r3790 r3983  
    2424class MyFCN : public FCNBase {
    2525public:
    26   MyFCN(double nu0,double dnu,TVector<r_8>& V,TVector<r_8>& EV);
     26  MyFCN(TVector<r_4>& FreqMHz,TVector<r_4>& V,TVector<r_4>& EV);
    2727  void InitFit(MnUserParameters& upar);
    28   double Func(int ix,const vector<double>& par) const;
     28  double Func(double fr,const vector<double>& par) const;
    2929  double operator()(const vector<double>& par) const;
    3030  double Up(void) const {return 1.;}
     31  void SetNu0(double nu0Mhz) {_nu0 = nu0Mhz;}
    3132  //protected:
     33  double _nu0;
    3234  mutable int _ndata;
    33   double _nu0, _dnu;
    34   TVector<r_8> _V;
    35   TVector<r_8> _EV;
     35  TVector<r_4> _F;
     36  TVector<r_4> _V;
     37  TVector<r_4> _EV;
    3638};
    3739} }  // namespace ROOT + Minuit2
    … …  
    4042int main(int narg,char *arg[])
    4143{
    42   double nu0 = 820.0, dnu = 0.5; // MHz
     44  int nslicemax = -1;
    4345  if(narg<=1) {
    44     cout<<"cmvgsmcfit cubefile.ppf"<<endl;
     46    cout<<"cmvgsmcfit cubefile.ppf [nslicemax]"<<endl;
    4547    return -1;
    4648  }
     49
     50  if(narg>2) nslicemax = atoi(arg[2]);
     51  cout<<"Fitting "<<arg[1]<<" , nslicemax="<<nslicemax<<endl;
     52
     53int rc = 0;
     54try {
    4755
    4856  POutPersist pos("cmvgsmcfit.ppf");
    … …  
    5058  TArray<r_4> Cube;
    5159  pis>>PPFNameTag("cube")>>Cube;
    52   TVector<r_8> V(Cube.SizeX()), EV(Cube.SizeX());
    53 
    54   const int nvar = 6;
    55   const char *vname[nvar] = {"x2r","ok","ph","th","A","N"};
     60  TVector<r_4> Freq;
     61  pis>>PPFNameTag("freq")>>Freq;
     62
     63  const int nvar = 8;
     64  const char *vname[nvar] = {"x2r","ok","ph","th","v0","vn","A","N"};
    5665  NTuple nt(nvar,vname);
    5766  double xnt[nvar]; for(int i=0;i<nvar;i++) xnt[i] = 0.;
    … …  
    5968  int nfit = 0;
    6069  for(int j=0;j<Cube.SizeY();j++) {
    61     cout<<"..."<<j<<endl;
     70    cout<<"...j = "<<j<<endl;
    6271  for(int k=0;k<Cube.SizeZ();k++) {
     72
     73    TVector<r_4> V(Cube.SizeX()), EV(Cube.SizeX());
    6374    for(int i=0;i<Cube.SizeX();i++) {
    6475      V(i) = Cube(i,j,k);
    … …  
    6677    }
    6778
    68     MyFCN myfcn(nu0,dnu,V,EV);
     79    MyFCN myfcn(Freq,V,EV);
     80    myfcn.SetNu0(Freq(0));
    6981    MnUserParameters upar;
    7082    myfcn.InitFit(upar);
    … …  
    89101    xnt[0] = xi2r; xnt[1] = (okfit)? 1.: 0.;
    90102    xnt[2] = j; xnt[3] = k;
    91     for(unsigned int i=0;i<par.size();i++) xnt[4+i] = par[i];
     103    xnt[4] = V(0); xnt[5] = V(V.Size()-1);
     104    for(unsigned int i=0;i<par.size();i++) xnt[6+i] = par[i];
    92105    nt.Fill(xnt);
    93106    //vector<double> err = ustate.Errors();
    94     for(int i=0;i<Cube.SizeX();i++) Cube(i,j,k) -= myfcn.Func(i,par);
     107    for(int i=0;i<Cube.SizeX();i++) Cube(i,j,k) -= myfcn.Func(Freq(i),par);
    95108
    96109    if(nfit<10) {
    97110      TVector<r_8> F(Cube.SizeX()), R(Cube.SizeX());
    98111      for(int i=0;i<Cube.SizeX();i++) {
    99         F(i) = myfcn.Func(i,par);
     112        F(i) = myfcn.Func(Freq(i),par);
    100113        R(i) = V(i) - F(i);
    101114      }
    … …  
    108121      pos<<PPFNameTag(str)<<R;
    109122    }
     123
    110124    nfit++;
    111     if(nfit>10) break;
     125
    112126  }
     127    if(nslicemax>0 && j>=nslicemax) break;
    113128  }
    114129
     130  cout<<"nfit = "<<nfit<<endl;
    115131  pos<<PPFNameTag("nt")<<nt;
    116132  pos<<PPFNameTag("cublis")<<Cube;
    117133
     134//-----------------------------------------------------------------
     135} catch (PException& exc) {
     136  cerr<<"cmvgsmcfit.cc catched PException"<<exc.Msg()<<endl;
     137  rc = 77;
     138} catch (std::exception& sex) {
     139  cerr << "cmvgsmcfit.cc std::exception :"
     140       << (string)typeid(sex).name() << "\n msg= " << sex.what() << endl;
     141  rc = 78;
     142} catch (...) {
     143  cerr << "cmvgsmcfit.cc catched unknown (...) exception  " << endl;
     144  rc = 79;
     145}
     146
     147return rc;
    118148}
    119149
    120150//--------------------------------------------------
    121 MyFCN::MyFCN(double nu0, double dnu,TVector<r_8>& V,TVector<r_8>& EV)
    122   : _V(V,false), _EV(EV,false)
    123 {
     151MyFCN::MyFCN(TVector<r_4>& FreqMHz,TVector<r_4>& V,TVector<r_4>& EV)
     152  : _F(FreqMHz,false), _V(V,false), _EV(EV,false)
     153{
     154  if(_F.Size()==0 || _F.Size()!=_V.Size() || _F.Size()!=_EV.Size())
     155    throw SzMismatchError("MyFCN::MyFCN_Error: incompatible F / V / EV sizes");
     156  _nu0 = _F(0);
    124157  _ndata =_V.Size();
    125   _nu0 = nu0;
    126   _dnu = dnu;
    127158}
    128159
    129160void MyFCN::InitFit(MnUserParameters& upar)
    130161{
    131   int n = _V.Size()-1;
     162  int n = _F.Size() - 1;
    132163  double a = _V(0);
    133   double ni = log(_V(n)/a) / log(1.+n*_dnu/_nu0);
     164  double ni = log(_V(n)/a) / log(_F(n)/_nu0);
    134165
    135166  upar.Add("A",a,a/50.);
    … …  
    137168}
    138169
    139 double MyFCN::Func(int ix,const vector<double>& par) const
    140 // f = A*(nu/nu0)^N   avec  nu = nu0 + i*dnu
    141 {
    142  double x = (_nu0 + ix*_dnu)/_nu0;
    143  double f = par[0]*pow(x,par[1]);
     170double MyFCN::Func(double fr,const vector<double>& par) const
     171// v = A*(nu/nu0)^N
     172{
     173 double f = par[0]*pow(fr/_nu0,par[1]);
    144174 return f;
    145175}
    … …  
    153183   if(e2<=0.) continue;
    154184   double v = _V(i);
    155    double f = Func(i,par);
     185   double f = Func(_F(i),par);
    156186   xi2 += (v-f)*(v-f)/e2;
    157187   _ndata++;
    … …  
    166196set j 0
    167197set k 0
     198zone 1 2
    168199n/plot v_${j}_${k}.val%n ! ! "nsta cpts"
    169200n/plot f_${j}_${k}.val%n ! ! "nsta cpts same red"
    170201n/plot r_${j}_${k}.val%n ! ! "nsta cpts"
    171202
    172 n/plot nt.ok%_nl ! ! "cpts"
    173 n/plot nt.x2r%_nl ! ! "cpts"
    174 
    175 n/plot nt.A%_nl ! ! "cpts"
    176 n/plot nt.N%_nl ! ! "cpts"
     203zone 1 2
     204n/plot nt.ok%_nl ! ! "crossmarker3"
     205n/plot nt.x2r%_nl ! ! "crossmarker3"
     206
     207zone 1 2
     208n/plot nt.v0%_nl ! ! "crossmarker3"
     209n/plot nt.vn%_nl ! ! "plusmarker3 same red"
     210n/plot nt.log10(v0/vn)%v0 ! ! "crossmarker3"
     211
     212
     213zone
     214n/plot nt.A%_nl ! ! "crossmarker3"
     215n/plot nt.N%_nl ! ! "crossmarker3"
     216
     217n/plot nt.A%v0 ! ! "crossmarker3"
     218n/plot nt.(A-v0)/v0%_nl ! ! "crossmarker3"
     219n/plot nt.N%A ! ! "crossmarker3"
    177220
    178221set n 0
Note: See TracChangeset for help on using the changeset viewer.