Index: trunk/SophyaLib/Samba/spherethetaphi.cc
===================================================================
--- trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 487)
+++ trunk/SophyaLib/Samba/spherethetaphi.cc	(revision 515)
@@ -4,4 +4,792 @@
 #include "piocmplx.h"
 #include <iostream.h>
+
+
+//***************************************************************
+//++ 
+// Class	SphereThetaPhi
+// 
+//
+// include      spherethetaphi.h nbmath.h 
+//
+//    Découpage de la sphère selon theta et phi, chaque 
+//    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du 
+//    beurre), chacune des m bandes de theta ainsi définies étant découpée par 
+//    des méridiens équirepartis, ce découpage étant fait de sorte que tous 
+//    les pixels aient la même surface et soient le plus carré possible.
+//    On commence par découper l'hémisphère de z positif en partant du pôle et 
+//    en allant vers l'équateur. Le premier pixel est la calotte polaire, 
+//    il est circulaire et centré sur theta=0.
+//--
+//++
+//
+// Links	Parents
+//
+//    SphericalMap 
+//-- 
+//++
+//
+// Links	Descendants
+//
+//     
+//-- 
+
+/* --Methode-- */
+//++
+// Titre	Constructeurs
+//--
+//++
+
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi()
+
+//--
+{
+  InitNul();
+}
+
+
+/* --Methode-- */
+
+//++
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(int m)
+
+//    Constructeur : m est le nombre de bandes en theta sur un hémisphère 
+//    (la calotte constituant la premiere bande).
+//    pet est le nombre de pixels (pétales) de la bande en contact avec la 
+//    calotte polaire. Pour l'instant pet est inopérant! 
+//--
+{
+  InitNul();
+  Pixelize(m);
+}
+
+template <class T>
+SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
+{
+  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
+  NTheta_= s.NTheta_;
+  NPix_  = s.NPix_;
+  NPhi_  = new int[NTheta_];
+  Theta_ = new double[NTheta_+1];
+  TNphi_ = new int[NTheta_+1];
+  for(int k = 0; k < NTheta_; k++) 
+    {
+      NPhi_[k] = s.NPhi_[k];
+      Theta_[k]= s.Theta_[k];
+      TNphi_[k]= s.TNphi_[k];
+    }
+  Theta_[NTheta_]= s.Theta_[NTheta_];
+  TNphi_[NTheta_]= s.TNphi_[NTheta_];
+  Omega_ = s.Omega_;
+  pixels_= s.pixels_;
+}
+
+//++
+// Titre	Destructeur
+//--
+//++
+template <class T>
+SphereThetaPhi<T>::~SphereThetaPhi()
+
+//--
+{
+  Clear();
+}
+
+//++
+// Titre	Méthodes
+//--
+template <class T>
+void SphereThetaPhi<T>::InitNul()
+//
+//    initialise à zéro les variables de classe pointeurs
+{
+  NTheta_= 0;
+  NPix_  = 0;
+  Theta_ = NULL;
+  NPhi_  = NULL;
+  TNphi_ = NULL;
+  pixels_.Reset();
+}
+
+/* --Methode-- */
+template <class T>
+void SphereThetaPhi<T>::Clear()
+
+{
+  if(Theta_)  delete[] Theta_;
+  if(NPhi_ )  delete[] NPhi_;
+  if(TNphi_)  delete[] TNphi_;
+  InitNul();
+}
+
+//++
+template <class T>
+void SphereThetaPhi<T>::Resize(int m)
+//   re-pixelize the sphere
+//--
+{
+  Clear();
+  Pixelize(m);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int SphereThetaPhi<T>::NbPixels() const
+
+//    Retourne le nombre de pixels du découpage
+//--
+{
+  return(NPix_);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+T& SphereThetaPhi<T>::PixVal(int k)
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if((k < 0) || (k >= NPix_)) 
+    {
+      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+      cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
+      THROW(rangeCheckErr);
+    }
+  return pixels_(k);
+}
+
+//++
+template <class T>
+T const& SphereThetaPhi<T>::PixVal(int k) const
+
+//    Retourne la valeur du contenu du pixel d'indice k
+//--
+{
+  if((k < 0) || (k >= NPix_)) 
+    {
+      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
+      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
+      throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
+    } 
+  return *(pixels_.Data()+k);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
+
+//  Retourne l'indice du pixel vers lequel pointe une direction définie par 
+//    ses coordonnées sphériques
+//--
+{
+  double dphi;
+  int i,j,k;
+  bool fgzn = false;
+
+  if((theta > Pi) || (theta < 0.)) return(-1);
+  if((phi < 0.) || (phi > DeuxPi)) return(-1);
+  if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
+  
+  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
+  // Theta_[kt] et Theta_[kt+1]
+  for( i=1; i< NTheta_; i++ ) 
+    if( theta < Theta_[i] ) break;
+  
+  dphi= DeuxPi/(double)NPhi_[i-1];
+  
+  if (fgzn)  k= NPix_-TNphi_[i]+(int)(phi/dphi);
+  else k= TNphi_[i-1]+(int)(phi/dphi); 
+  return(k);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const
+
+//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
+//--
+{
+  int i; 
+  bool fgzn = false;
+  
+  if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
+  if( k >= NPix_/2)  {fgzn = true; k = NPix_-1-k;}
+
+  // recupère l'indice i de la tranche qui contient le pixel k
+  for( i=0; i< NTheta_; i++ ) 
+    if( k < TNphi_[i+1] ) break;
+
+  // angle theta
+  theta= 0.5*(Theta_[i]+Theta_[i+1]);
+  if (fgzn) theta= Pi-theta;
+  
+  // angle phi
+  k -= TNphi_[i];
+  phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);
+  if (fgzn) phi= DeuxPi-phi;
+}
+
+//++
+template <class T>
+double SphereThetaPhi<T>::PixSolAngle(int dummy) const
+
+//    Pixel Solid angle  (steradians)
+//    All the pixels have the same solid angle. The dummy argument is
+//    for compatibility with eventual pixelizations which would not 
+//    fulfil this requirement.
+//--
+{
+  return Omega_;
+}
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax) 
+
+//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
+//--
+  {
+  int j;
+  double dphi;
+  bool fgzn= false;
+  
+  if((k < 0) || (k >= NPix_)) {
+    tetMin= -99999.; 
+    phiMin= -99999.;  
+    tetMax= -99999.; 
+    phiMax= -99999.;  
+    return;
+  }
+  
+  // si on se trouve dans l'hémisphère Sud
+  if(k >= NPix_/2) { 
+    fgzn= true; 
+    k= NPix_-1-k; 
+  }
+  
+  // on recupere l'indice i de la tranche qui contient le pixel k
+  int i;
+  for( i=0; i< NTheta_; i++ ) 
+    if(k < TNphi_[i+1]) break;
+  
+  // valeurs limites de theta dans l'hemisphere Nord
+  tetMin= Theta_[i];
+  tetMax= Theta_[i+1];
+  // valeurs limites de theta dans l'hemisphere Sud
+  if (fgzn) {
+    tetMin= Pi-Theta_[i+1];
+    tetMax= Pi-Theta_[i];
+  }
+  
+  // pixel correspondant dans l'hemisphere Nord
+  if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
+  
+  // indice j de discretisation ( phi= j*dphi )
+  j= k-TNphi_[i];
+  dphi= DeuxPi/(double)NPhi_[i];
+  
+  // valeurs limites de phi 
+  phiMin= dphi*(double)(j);
+  phiMax= dphi*(double)(j+1);
+  return;
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int SphereThetaPhi<T>::NbThetaSlices() const
+
+//    Retourne le nombre de tranches en theta sur la sphere
+//--
+{
+  int nbslices;
+  nbslices= 2*NTheta_;
+  return(nbslices);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int SphereThetaPhi<T>::NPhi(int kt) const
+
+//    Retourne le nombre de pixels en phi de la tranche kt
+//-- 
+{
+  int nbpix;  
+  // verification
+  if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
+  
+  // si on se trouve dans l'hemisphere Sud
+  if(kt >= NTheta_) {
+    kt= 2*NTheta_-1-kt; 
+  }
+  
+  // nombre de pixels
+  nbpix= NPhi_[kt];
+  return(nbpix);
+}
+
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax) 
+
+//    Retourne les valeurs de theta limitant la tranche kt
+//--
+{
+  bool fgzn= false;
+  // verification
+  if( (kt< 0) || (kt>= 2*NTheta_) ) {
+    tetMin= -99999.; 
+    tetMax= -99999.; 
+    return;
+  }
+
+  // si on se trouve dans l'hemisphere Sud
+  if( kt >= NTheta_ ) {
+    fgzn= true;
+    kt= 2*NTheta_-1-kt; 
+  }
+  
+  // valeurs limites de theta dans l'hemisphere Nord
+  tetMin= Theta_[kt];
+  tetMax= Theta_[kt+1];
+  // valeurs limites de theta dans l'hemisphere Sud
+  if (fgzn) {
+    tetMin= Pi-Theta_[kt+1];
+    tetMax= Pi-Theta_[kt];
+  }  
+}
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax) 
+
+//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
+//--
+{
+  // verification
+  if((kt < 0) || (kt >= 2*NTheta_)) {
+    phiMin= -99999.; 
+    phiMax= -99999.; 
+    return;
+  }
+  
+  // si on se trouve dans l'hemisphere Sud
+  if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
+  
+  // verifie si la tranche kt contient au moins jp pixels
+  if( (jp< 0) || (jp >= NPhi_[kt]) ) {
+    phiMin= -88888.; 
+    phiMax= -88888.; 
+    return;
+  }
+  
+  double dphi= DeuxPi/(double)NPhi_[kt];
+  phiMin= dphi*(double)(jp);
+  phiMax= dphi*(double)(jp+1);
+  return;
+}
+
+/* --Methode-- */
+//++
+template <class T>
+int SphereThetaPhi<T>::Index(int kt,int jp) const
+
+//    Retourne l'indice du pixel d'indice jp dans la tranche kt 
+//--
+{
+  int k;
+  bool fgzn= false;
+  
+  // si on se trouve dans l'hemisphere Sud
+  if(kt >= NTheta_) {
+    fgzn= true; 
+    kt= 2*NTheta_-1-kt;
+  }
+  
+  // si la tranche kt contient au moins i pixels
+  if( (jp>=0) && (jp<NPhi_[kt]) ) 
+    {
+      // dans l'hemisphere Sud
+      if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
+      // dans l'hemisphere Nord
+      else k= TNphi_[kt]+jp;
+    }
+  else
+    {
+      k= 9999;
+      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
+    }
+  return(k);
+}
+
+/* --Methode-- */
+//++
+template <class T>
+void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp) 
+
+//    Retourne les indices kt et jp du pixel d'indice k
+//--
+{
+  bool fgzn= false;  
+  // si on se trouve dans l'hemisphere Sud
+  if(k >= NPix_/2) 
+    { 
+      fgzn= true; 
+      k= NPix_-1-k; 
+    }
+  
+  // on recupere l'indice kt de la tranche qui contient le pixel k
+  int i;
+  for(i = 0; i < NTheta_; i++) 
+    if(k < TNphi_[i+1]) break;
+  
+  // indice  kt de tranche
+  if (fgzn) kt= 2*NTheta_-1-i;
+  else kt= i;
+  
+  // indice jp de pixel
+  if (fgzn) jp= TNphi_[i+1]-k-1;
+  else jp= k-TNphi_[i];  
+}
+//++
+template <class T>
+void  SphereThetaPhi<T>::Pixelize(int m) 
+
+//    effectue le découpage en pixels (m et pet ont la même signification 
+//    que pour le constructeur)
+//
+//    Chaque bande de theta sera découpée en partant de phi=0 ...
+//    L'autre hémisphère est parcourue dans le même sens en phi et de
+//    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
+//    proche de l'équateur a z>0 est celui de plus petit phi de la bande la 
+//    plus proche de l'equateur a z<0).
+//--
+{
+  int ntotpix,i,j;
+  
+  // Decodage et controle des arguments d'appel :
+  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
+  if (m < 2) m = 2;
+  if (m > 16384) m = 16384;
+  
+  // On memorise les arguments d'appel
+  NTheta_ = m;  
+  
+  // On commence par decouper l'hemisphere z>0.
+  // Creation des vecteurs contenant :
+  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
+  Theta_= new double[m+1];
+  
+  // Le nombre de pixels en phi de chacune des bandes en theta
+  NPhi_ = new int[m];
+  
+  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
+  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
+  TNphi_= new int[m+1];
+  
+  // Calcul du nombre total de pixels dans chaque bande pour optimiser
+  // le rapport largeur/hauteur des pixels 
+  
+  //calotte polaire
+  TNphi_[0]= 0;
+  NPhi_[0] = 1;
+  
+  //bandes jusqu'a l'equateur
+  for(j = 1; j < m; j++) 
+    {
+      TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
+      NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
+    }
+  
+  // Nombre total de pixels sur l'hemisphere
+  ntotpix  = TNphi_[m-1]+NPhi_[m-1];
+  TNphi_[m]= ntotpix;
+  // et sur la sphere entiere
+  NPix_= 2*ntotpix;
+  
+  // Creation et initialisation du vecteur des contenus des pixels 
+  pixels_.ReSize(NPix_);
+  pixels_.Reset();
+
+  // Determination des limites des bandes en theta :
+  // omeg est l'angle solide couvert par chaque pixel,
+  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
+  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide 
+  //de la
+  // calotte allant du pole a la limite haute de la bande kt vaut
+  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
+  
+  double omeg2pi= 1./(double)ntotpix;
+  Omega_ = omeg2pi*DeuxPi;
+  
+  for(j=0; j <= m; j++) 
+    {
+      Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
+    }
+}
+
+//++
+template <class T>
+void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const 
+
+//    Retourne, pour la tranche en theta d'indice 'index' le theta 
+//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
+//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
+//    correspondantes
+//--
+
+{
+  cout << "entree GetThetaSlice, couche no " << index << endl;
+
+  if(index < 0 || index > NbThetaSlices()) 
+    {
+      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
+      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
+      THROW(rangeCheckErr);
+    }
+
+  int iring= Index(index,0);
+  int bid  = this->NPhi(index);
+  int lring  = bid;
+  cout << " iring= " << iring << " lring= " << lring << endl;
+
+  phi.ReSize(lring);
+  value.ReSize(lring);
+  double Te= 0.;
+  double Fi= 0.;
+  for(int kk = 0; kk < lring; kk++) 
+    {
+      PixThetaPhi(kk+iring,Te,Fi);
+      phi(kk)= Fi;
+      value(kk)= PixVal(kk+iring);
+    }
+  theta= Te;
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmNPhi(int* array, int m)
+  //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta 
+  //--
+{
+  NPhi_= new int[m];
+  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmTNphi(int* array, int m)
+  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes 
+  //--
+{
+  TNphi_= new int[m];
+  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setmTheta(double* array, int m)
+  //remplit  le tableau contenant les valeurs limites de theta
+  //--
+{
+  Theta_= new double[m];
+  for(int k = 0; k < m; k++) Theta_[k]= array[k];
+}
+
+template <class T>
+void SphereThetaPhi<T>::setDataBlock(T* data, int m)
+  // remplit  le vecteur des contenus des pixels 
+{
+  pixels_.FillFrom(m,data);
+} 
+
+template <class T>
+void SphereThetaPhi<T>::print(ostream& os) const
+{
+  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
+  //
+  os << " NTheta_= " << NTheta_ << endl;
+  os << " NPix_    = " << NPix_   << endl;
+  os << " Omega_  =  " << Omega_   << endl;
+
+  os << " contenu de NPhi_ : ";
+  int i;
+  for(i=0; i < NTheta_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << NPhi_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de Theta_ : ";
+  for(i=0; i < NTheta_+1; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << Theta_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de TNphi_ : ";
+  for(i=0; i < NTheta_+1; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os << TNphi_[i] <<", ";
+    }
+  os << endl;
+
+  os << " contenu de pixels : ";
+  for(i=0; i < NPix_; i++)
+    {
+      if(i%5 == 0) os << endl;
+      os <<  pixels_(i) <<", ";
+    }
+  os << endl;
+}
+
+///////////////////////////////////////////////////////////
+// --------------------------------------------------------
+//   Les objets delegues pour la gestion de persistance 
+// --------------------------------------------------------
+//////////////////////////////////////////////////////////
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
+{
+  dobj= new SphereThetaPhi<T>;
+  ownobj= true;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
+{
+  dobj= new SphereThetaPhi<T>;
+  ownobj= true;
+  Read(filename);
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
+{
+  dobj= new SphereThetaPhi<T>(obj);
+  ownobj= true;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
+{
+  dobj= obj;
+  ownobj= false;
+}
+
+template <class T>
+FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
+{
+  if (ownobj && dobj) delete dobj;
+}
+
+template <class T>
+AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
+{
+  return(dobj);
+}
+
+template <class T>
+void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
+{
+  cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      dobj= new SphereThetaPhi<T>;
+    }
+
+  // Pour savoir s'il y avait un DVList Info associe
+  char strg[256];
+  is.GetLine(strg, 255);
+  bool hadinfo= false;
+  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
+  if(hadinfo) 
+    {    // Lecture eventuelle du DVList Info
+      is >> dobj->Info(); 
+    }
+
+  int mNTheta;
+  is.GetI4(mNTheta);  
+  dobj->setSizeIndex(mNTheta);
+
+  int mNPix;
+  is.GetI4(mNPix);
+  dobj->setNbPixels(mNPix);
+
+  double mOmeg;
+  is.GetR8(mOmeg);
+  dobj->setPixSolAngle(mOmeg);
+
+  int* mNphi= new int[mNTheta];
+  is.GetI4s(mNphi, mNTheta);
+  dobj->setmNPhi(mNphi, mNTheta);
+  delete [] mNphi;
+
+  int* mTNphi= new int[mNTheta+1];
+  is.GetI4s(mTNphi, mNTheta+1);
+  dobj->setmTNphi(mTNphi, mNTheta+1);
+  delete [] mTNphi;
+
+  double* mTheta= new double[mNTheta+1];
+  is.GetR8s(mTheta, mNTheta+1);
+  dobj->setmTheta(mTheta, mNTheta+1);
+  delete [] mTheta;
+
+  T* mPixels= new T[mNPix];
+  PIOSReadArray(is, mPixels, mNPix);
+  dobj->setDataBlock(mPixels, mNPix);
+  delete [] mPixels;
+}
+
+template <class T>
+void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
+{
+  cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
+
+  if(dobj == NULL) 
+    {
+      cout << " WriteSelf:: dobj= null " << endl;
+      return;
+    }
+
+  char strg[256];
+  int mNTheta= dobj->SizeIndex();
+  int mNPix  = dobj->NbPixels();
+ 
+  if(dobj->ptrInfo()) 
+    {
+      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
+      os.PutLine(strg);
+      os << dobj->Info();
+    }
+  else 
+    { 
+      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
+      os.PutLine(strg);  
+    }
+
+  os.PutI4(mNTheta);  
+  os.PutI4(mNPix);
+  os.PutR8(dobj->PixSolAngle(0));
+  os.PutI4s(dobj->getmNPhi() , mNTheta);
+  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
+  os.PutR8s(dobj->getmTheta(), mNTheta+1);
+  //os.Put((dobj->getDataBlock())->Data(), mNPix);
+  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
+}
 
 #ifdef __CXX_PRAGMA_TEMPLATES__
@@ -25,788 +813,2 @@
 template class FIO_SphereThetaPhi< complex<double> >;
 #endif
-
-//***************************************************************
-//++ 
-// Class	SphereThetaPhi
-// 
-//
-// include      spherethetaphi.h nbmath.h 
-//
-//    Découpage de la sphère selon theta et phi, chaque 
-//    hémisphère étant découpé en (m-1) parallèles (l'équateur compte pour du 
-//    beurre), chacune des m bandes de theta ainsi définies étant découpée par 
-//    des méridiens équirepartis, ce découpage étant fait de sorte que tous 
-//    les pixels aient la même surface et soient le plus carré possible.
-//    On commence par découper l'hémisphère de z positif en partant du pôle et 
-//    en allant vers l'équateur. Le premier pixel est la calotte polaire, 
-//    il est circulaire et centré sur theta=0.
-//--
-//++
-//
-// Links	Parents
-//
-//    SphericalMap 
-//-- 
-//++
-//
-// Links	Descendants
-//
-//     
-//-- 
-
-/* --Methode-- */
-//++
-// Titre	Constructeurs
-//--
-//++
-
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi()
-
-//--
-{
-  InitNul();
-}
-
-
-/* --Methode-- */
-
-//++
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi(int m)
-
-//    Constructeur : m est le nombre de bandes en theta sur un hémisphère 
-//    (la calotte constituant la premiere bande).
-//    pet est le nombre de pixels (pétales) de la bande en contact avec la 
-//    calotte polaire. Pour l'instant pet est inopérant! 
-//--
-{
-  InitNul();
-  Pixelize(m);
-}
-
-template <class T>
-SphereThetaPhi<T>::SphereThetaPhi(const SphereThetaPhi<T>& s)
-{
-  if(s.mInfo_) mInfo_= new DVList(*s.mInfo_);
-  NTheta_= s.NTheta_;
-  NPix_  = s.NPix_;
-  NPhi_  = new int[NTheta_];
-  Theta_ = new double[NTheta_+1];
-  TNphi_ = new int[NTheta_+1];
-  for(int k = 0; k < NTheta_; k++) 
-    {
-      NPhi_[k] = s.NPhi_[k];
-      Theta_[k]= s.Theta_[k];
-      TNphi_[k]= s.TNphi_[k];
-    }
-  Theta_[NTheta_]= s.Theta_[NTheta_];
-  TNphi_[NTheta_]= s.TNphi_[NTheta_];
-  Omega_ = s.Omega_;
-  pixels_= s.pixels_;
-}
-
-//++
-// Titre	Destructeur
-//--
-//++
-template <class T>
-SphereThetaPhi<T>::~SphereThetaPhi()
-
-//--
-{
-  Clear();
-}
-
-//++
-// Titre	Méthodes
-//--
-template <class T>
-void SphereThetaPhi<T>::InitNul()
-//
-//    initialise à zéro les variables de classe pointeurs
-{
-  NTheta_= 0;
-  NPix_  = 0;
-  Theta_ = NULL;
-  NPhi_  = NULL;
-  TNphi_ = NULL;
-  pixels_.Reset();
-}
-
-/* --Methode-- */
-template <class T>
-void SphereThetaPhi<T>::Clear()
-
-{
-  if(Theta_)  delete[] Theta_;
-  if(NPhi_ )  delete[] NPhi_;
-  if(TNphi_)  delete[] TNphi_;
-  InitNul();
-}
-
-//++
-template <class T>
-void SphereThetaPhi<T>::Resize(int m)
-//   re-pixelize the sphere
-//--
-{
-  Clear();
-  Pixelize(m);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int SphereThetaPhi<T>::NbPixels() const
-
-//    Retourne le nombre de pixels du découpage
-//--
-{
-  return(NPix_);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-T& SphereThetaPhi<T>::PixVal(int k)
-
-//    Retourne la valeur du contenu du pixel d'indice k
-//--
-{
-  if((k < 0) || (k >= NPix_)) 
-    {
-      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
-      cout << " SphereThetaPhi::PIxVal : exceptions a mettre en place" <<endl;
-      THROW(rangeCheckErr);
-    }
-  return pixels_(k);
-}
-
-//++
-template <class T>
-T const& SphereThetaPhi<T>::PixVal(int k) const
-
-//    Retourne la valeur du contenu du pixel d'indice k
-//--
-{
-  if((k < 0) || (k >= NPix_)) 
-    {
-      cout << " SphereThetaPhi::PIxVal : exceptions  a mettre en place" <<endl;
-      //THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range"));
-      throw "SphereThetaPhi::PIxVal Pixel index out of range"; 
-    } 
-  return *(pixels_.Data()+k);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int SphereThetaPhi<T>::PixIndexSph(double theta, double phi) const
-
-//  Retourne l'indice du pixel vers lequel pointe une direction définie par 
-//    ses coordonnées sphériques
-//--
-{
-  double dphi;
-  int i,j,k;
-  bool fgzn = false;
-
-  if((theta > Pi) || (theta < 0.)) return(-1);
-  if((phi < 0.) || (phi > DeuxPi)) return(-1);
-  if(theta > Pi*0.5) {fgzn = true; theta = Pi-theta;}
-  
-  // La bande d'indice kt est limitée par les valeurs de theta contenues dans
-  // Theta_[kt] et Theta_[kt+1]
-  for( i=1; i< NTheta_; i++ ) 
-    if( theta < Theta_[i] ) break;
-  
-  dphi= DeuxPi/(double)NPhi_[i-1];
-  
-  if (fgzn)  k= NPix_-TNphi_[i]+(int)(phi/dphi);
-  else k= TNphi_[i-1]+(int)(phi/dphi); 
-  return(k);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::PixThetaPhi(int k,double& theta,double& phi) const
-
-//    Retourne les coordonnées (theta,phi) du milieu du pixel d'indice k
-//--
-{
-  int i; 
-  bool fgzn = false;
-  
-  if((k < 0) || (k >= NPix_)) {theta = -99999.; phi = -99999.; return; }
-  if( k >= NPix_/2)  {fgzn = true; k = NPix_-1-k;}
-
-  // recupère l'indice i de la tranche qui contient le pixel k
-  for( i=0; i< NTheta_; i++ ) 
-    if( k < TNphi_[i+1] ) break;
-
-  // angle theta
-  theta= 0.5*(Theta_[i]+Theta_[i+1]);
-  if (fgzn) theta= Pi-theta;
-  
-  // angle phi
-  k -= TNphi_[i];
-  phi= DeuxPi/(double)NPhi_[i]*(double)(k+.5);
-  if (fgzn) phi= DeuxPi-phi;
-}
-
-//++
-template <class T>
-double SphereThetaPhi<T>::PixSolAngle(int dummy) const
-
-//    Pixel Solid angle  (steradians)
-//    All the pixels have the same solid angle. The dummy argument is
-//    for compatibility with eventual pixelizations which would not 
-//    fulfil this requirement.
-//--
-{
-  return Omega_;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Limits(int k,double& tetMin,double& tetMax,double& phiMin,double& phiMax) 
-
-//    Retourne les valeurs de theta et phi limitant le pixel d'indice k
-//--
-  {
-  int j;
-  double dphi;
-  bool fgzn= false;
-  
-  if((k < 0) || (k >= NPix_)) {
-    tetMin= -99999.; 
-    phiMin= -99999.;  
-    tetMax= -99999.; 
-    phiMax= -99999.;  
-    return;
-  }
-  
-  // si on se trouve dans l'hémisphère Sud
-  if(k >= NPix_/2) { 
-    fgzn= true; 
-    k= NPix_-1-k; 
-  }
-  
-  // on recupere l'indice i de la tranche qui contient le pixel k
-  int i;
-  for( i=0; i< NTheta_; i++ ) 
-    if(k < TNphi_[i+1]) break;
-  
-  // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= Theta_[i];
-  tetMax= Theta_[i+1];
-  // valeurs limites de theta dans l'hemisphere Sud
-  if (fgzn) {
-    tetMin= Pi-Theta_[i+1];
-    tetMax= Pi-Theta_[i];
-  }
-  
-  // pixel correspondant dans l'hemisphere Nord
-  if (fgzn) k= TNphi_[i+1]-k+TNphi_[i]-1;
-  
-  // indice j de discretisation ( phi= j*dphi )
-  j= k-TNphi_[i];
-  dphi= DeuxPi/(double)NPhi_[i];
-  
-  // valeurs limites de phi 
-  phiMin= dphi*(double)(j);
-  phiMax= dphi*(double)(j+1);
-  return;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int SphereThetaPhi<T>::NbThetaSlices() const
-
-//    Retourne le nombre de tranches en theta sur la sphere
-//--
-{
-  int nbslices;
-  nbslices= 2*NTheta_;
-  return(nbslices);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int SphereThetaPhi<T>::NPhi(int kt) const
-
-//    Retourne le nombre de pixels en phi de la tranche kt
-//-- 
-{
-  int nbpix;  
-  // verification
-  if((kt < 0) || (kt >= 2*NTheta_)) return(-1);
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) {
-    kt= 2*NTheta_-1-kt; 
-  }
-  
-  // nombre de pixels
-  nbpix= NPhi_[kt];
-  return(nbpix);
-}
-
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Theta(int kt,double& tetMin,double& tetMax) 
-
-//    Retourne les valeurs de theta limitant la tranche kt
-//--
-{
-  bool fgzn= false;
-  // verification
-  if( (kt< 0) || (kt>= 2*NTheta_) ) {
-    tetMin= -99999.; 
-    tetMax= -99999.; 
-    return;
-  }
-
-  // si on se trouve dans l'hemisphere Sud
-  if( kt >= NTheta_ ) {
-    fgzn= true;
-    kt= 2*NTheta_-1-kt; 
-  }
-  
-  // valeurs limites de theta dans l'hemisphere Nord
-  tetMin= Theta_[kt];
-  tetMax= Theta_[kt+1];
-  // valeurs limites de theta dans l'hemisphere Sud
-  if (fgzn) {
-    tetMin= Pi-Theta_[kt+1];
-    tetMax= Pi-Theta_[kt];
-  }  
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::Phi(int kt,int jp,double& phiMin,double& phiMax) 
-
-//    Retourne les valeurs de phi limitant le pixel jp de la tranche kt
-//--
-{
-  // verification
-  if((kt < 0) || (kt >= 2*NTheta_)) {
-    phiMin= -99999.; 
-    phiMax= -99999.; 
-    return;
-  }
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) kt= 2*NTheta_-1-kt;
-  
-  // verifie si la tranche kt contient au moins jp pixels
-  if( (jp< 0) || (jp >= NPhi_[kt]) ) {
-    phiMin= -88888.; 
-    phiMax= -88888.; 
-    return;
-  }
-  
-  double dphi= DeuxPi/(double)NPhi_[kt];
-  phiMin= dphi*(double)(jp);
-  phiMax= dphi*(double)(jp+1);
-  return;
-}
-
-/* --Methode-- */
-//++
-template <class T>
-int SphereThetaPhi<T>::Index(int kt,int jp) const
-
-//    Retourne l'indice du pixel d'indice jp dans la tranche kt 
-//--
-{
-  int k;
-  bool fgzn= false;
-  
-  // si on se trouve dans l'hemisphere Sud
-  if(kt >= NTheta_) {
-    fgzn= true; 
-    kt= 2*NTheta_-1-kt;
-  }
-  
-  // si la tranche kt contient au moins i pixels
-  if( (jp>=0) && (jp<NPhi_[kt]) ) 
-    {
-      // dans l'hemisphere Sud
-      if (fgzn) k= NPix_-TNphi_[kt+1]+jp;
-      // dans l'hemisphere Nord
-      else k= TNphi_[kt]+jp;
-    }
-  else
-    {
-      k= 9999;
-      printf("\n la tranche %d ne contient pas un pixel de rang %d",kt,jp);
-    }
-  return(k);
-}
-
-/* --Methode-- */
-//++
-template <class T>
-void SphereThetaPhi<T>::ThetaPhiIndex(int k,int& kt,int& jp) 
-
-//    Retourne les indices kt et jp du pixel d'indice k
-//--
-{
-  bool fgzn= false;  
-  // si on se trouve dans l'hemisphere Sud
-  if(k >= NPix_/2) 
-    { 
-      fgzn= true; 
-      k= NPix_-1-k; 
-    }
-  
-  // on recupere l'indice kt de la tranche qui contient le pixel k
-  int i;
-  for(i = 0; i < NTheta_; i++) 
-    if(k < TNphi_[i+1]) break;
-  
-  // indice  kt de tranche
-  if (fgzn) kt= 2*NTheta_-1-i;
-  else kt= i;
-  
-  // indice jp de pixel
-  if (fgzn) jp= TNphi_[i+1]-k-1;
-  else jp= k-TNphi_[i];  
-}
-//++
-template <class T>
-void  SphereThetaPhi<T>::Pixelize(int m) 
-
-//    effectue le découpage en pixels (m et pet ont la même signification 
-//    que pour le constructeur)
-//
-//    Chaque bande de theta sera découpée en partant de phi=0 ...
-//    L'autre hémisphère est parcourue dans le même sens en phi et de
-//    l'équateur vers le pôle (le pixel qui suit le dernier de la bande la plus
-//    proche de l'équateur a z>0 est celui de plus petit phi de la bande la 
-//    plus proche de l'equateur a z<0).
-//--
-{
-  int ntotpix,i,j;
-  
-  // Decodage et controle des arguments d'appel :
-  // au moins 2 et au plus 16384 bandes d'un hemisphere en theta
-  if (m < 2) m = 2;
-  if (m > 16384) m = 16384;
-  
-  // On memorise les arguments d'appel
-  NTheta_ = m;  
-  
-  // On commence par decouper l'hemisphere z>0.
-  // Creation des vecteurs contenant :
-  // Les valeurs limites de theta (une valeur de plus que le nombre de bandes...)
-  Theta_= new double[m+1];
-  
-  // Le nombre de pixels en phi de chacune des bandes en theta
-  NPhi_ = new int[m];
-  
-  // Le nombre/Deuxpi total des pixels contenus dans les bandes de z superieur a une
-  // bande donnee (mTPphi[m] contient le nombre de pixels total de l'hemisphere)
-  TNphi_= new int[m+1];
-  
-  // Calcul du nombre total de pixels dans chaque bande pour optimiser
-  // le rapport largeur/hauteur des pixels 
-  
-  //calotte polaire
-  TNphi_[0]= 0;
-  NPhi_[0] = 1;
-  
-  //bandes jusqu'a l'equateur
-  for(j = 1; j < m; j++) 
-    {
-      TNphi_[j]= TNphi_[j-1]+NPhi_[j-1];
-      NPhi_[j] = (int)(.5+4.*(double)(m-.5)*sin(Pi*(double)j/(double)(2.*m-1.))) ;
-    }
-  
-  // Nombre total de pixels sur l'hemisphere
-  ntotpix  = TNphi_[m-1]+NPhi_[m-1];
-  TNphi_[m]= ntotpix;
-  // et sur la sphere entiere
-  NPix_= 2*ntotpix;
-  
-  // Creation et initialisation du vecteur des contenus des pixels 
-  pixels_.ReSize(NPix_);
-  pixels_.Reset();
-
-  // Determination des limites des bandes en theta :
-  // omeg est l'angle solide couvert par chaque pixel,
-  // une bande donnee kt couvre un angle solide NPhi_[kt]*omeg
-  // egal a 2* Pi*(cos Theta_[kt]-cos Theta_[kt+1]). De meme, l'angle solide 
-  //de la
-  // calotte allant du pole a la limite haute de la bande kt vaut
-  // 2* Pi*(1.-cos Theta_[kt+1])= TNphi_[kt]*omeg...
-  
-  double omeg2pi= 1./(double)ntotpix;
-  Omega_ = omeg2pi*DeuxPi;
-  
-  for(j=0; j <= m; j++) 
-    {
-      Theta_[j]= acos(1.-(double)TNphi_[j]*omeg2pi);
-    }
-}
-
-//++
-template <class T>
-void SphereThetaPhi<T>::GetThetaSlice(int index,double& theta, TVector<double>& phi, TVector<T>& value) const 
-
-//    Retourne, pour la tranche en theta d'indice 'index' le theta 
-//    correspondant, un vecteur (Peida) contenant les phi des pixels de 
-//    la tranche, un vecteur (Peida) contenant les valeurs de pixel 
-//    correspondantes
-//--
-
-{
-  cout << "entree GetThetaSlice, couche no " << index << endl;
-
-  if(index < 0 || index > NbThetaSlices()) 
-    {
-      // THROW(out_of_range("SphereThetaPhi::PIxVal Pixel index out of range")); 
-      cout << " SphereThetaPhi::GetThetaSlice : exceptions  a mettre en place" <<endl;
-      THROW(rangeCheckErr);
-    }
-
-  int iring= Index(index,0);
-  int bid  = this->NPhi(index);
-  int lring  = bid;
-  cout << " iring= " << iring << " lring= " << lring << endl;
-
-  phi.ReSize(lring);
-  value.ReSize(lring);
-  double Te= 0.;
-  double Fi= 0.;
-  for(int kk = 0; kk < lring; kk++) 
-    {
-      PixThetaPhi(kk+iring,Te,Fi);
-      phi(kk)= Fi;
-      value(kk)= PixVal(kk+iring);
-    }
-  theta= Te;
-}
-
-template <class T>
-void SphereThetaPhi<T>::setmNPhi(int* array, int m)
-  //remplit le tableau contenant le nombre de pixels en phi de chacune des bandes en theta 
-  //--
-{
-  NPhi_= new int[m];
-  for(int k = 0; k < m; k++) NPhi_[k]= array[k];
-}
-
-template <class T>
-void SphereThetaPhi<T>::setmTNphi(int* array, int m)
-  //remplit  le tableau contenant le nombre/Deuxpi total des pixels contenus dans les bandes 
-  //--
-{
-  TNphi_= new int[m];
-  for(int k = 0; k < m; k++) TNphi_[k]= array[k];
-}
-
-template <class T>
-void SphereThetaPhi<T>::setmTheta(double* array, int m)
-  //remplit  le tableau contenant les valeurs limites de theta
-  //--
-{
-  Theta_= new double[m];
-  for(int k = 0; k < m; k++) Theta_[k]= array[k];
-}
-
-template <class T>
-void SphereThetaPhi<T>::setDataBlock(T* data, int m)
-  // remplit  le vecteur des contenus des pixels 
-{
-  pixels_.FillFrom(m,data);
-} 
-
-template <class T>
-void SphereThetaPhi<T>::print(ostream& os) const
-{
-  if(mInfo_) os << "  DVList Info= " << *mInfo_ << endl;
-  //
-  os << " NTheta_= " << NTheta_ << endl;
-  os << " NPix_    = " << NPix_   << endl;
-  os << " Omega_  =  " << Omega_   << endl;
-
-  os << " contenu de NPhi_ : ";
-  for(int i=0; i < NTheta_; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << NPhi_[i] <<", ";
-    }
-  os << endl;
-
-  os << " contenu de Theta_ : ";
-  for(int i=0; i < NTheta_+1; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << Theta_[i] <<", ";
-    }
-  os << endl;
-
-  os << " contenu de TNphi_ : ";
-  for(int i=0; i < NTheta_+1; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os << TNphi_[i] <<", ";
-    }
-  os << endl;
-
-  os << " contenu de pixels : ";
-  for(int i=0; i < NPix_; i++)
-    {
-      if(i%5 == 0) os << endl;
-      os <<  pixels_(i) <<", ";
-    }
-  os << endl;
-}
-
-///////////////////////////////////////////////////////////
-// --------------------------------------------------------
-//   Les objets delegues pour la gestion de persistance 
-// --------------------------------------------------------
-//////////////////////////////////////////////////////////
-
-template <class T>
-FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi()
-{
-  dobj= new SphereThetaPhi<T>;
-  ownobj= true;
-}
-
-template <class T>
-FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(string const& filename)
-{
-  dobj= new SphereThetaPhi<T>;
-  ownobj= true;
-  Read(filename);
-}
-
-template <class T>
-FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(const SphereThetaPhi<T>& obj)
-{
-  dobj= new SphereThetaPhi<T>(obj);
-  ownobj= true;
-}
-
-template <class T>
-FIO_SphereThetaPhi<T>::FIO_SphereThetaPhi(SphereThetaPhi<T>* obj)
-{
-  dobj= obj;
-  ownobj= false;
-}
-
-template <class T>
-FIO_SphereThetaPhi<T>::~FIO_SphereThetaPhi()
-{
-  if (ownobj && dobj) delete dobj;
-}
-
-template <class T>
-AnyDataObj* FIO_SphereThetaPhi<T>::DataObj()
-{
-  return(dobj);
-}
-
-template <class T>
-void FIO_SphereThetaPhi<T>::ReadSelf(PInPersist& is)
-{
-  cout << " FIO_SphereThetaPhi:: ReadSelf " << endl;
-
-  if(dobj == NULL) 
-    {
-      dobj= new SphereThetaPhi<T>;
-    }
-
-  // Pour savoir s'il y avait un DVList Info associe
-  char strg[256];
-  is.GetLine(strg, 255);
-  bool hadinfo= false;
-  if(strncmp(strg+strlen(strg)-7, "HasInfo", 7) == 0)  hadinfo= true;
-  if(hadinfo) 
-    {    // Lecture eventuelle du DVList Info
-      is >> dobj->Info(); 
-    }
-
-  int mNTheta;
-  is.GetI4(mNTheta);  
-  dobj->setSizeIndex(mNTheta);
-
-  int mNPix;
-  is.GetI4(mNPix);
-  dobj->setNbPixels(mNPix);
-
-  double mOmeg;
-  is.GetR8(mOmeg);
-  dobj->setPixSolAngle(mOmeg);
-
-  int* mNphi= new int[mNTheta];
-  is.GetI4s(mNphi, mNTheta);
-  dobj->setmNPhi(mNphi, mNTheta);
-  delete [] mNphi;
-
-  int* mTNphi= new int[mNTheta+1];
-  is.GetI4s(mTNphi, mNTheta+1);
-  dobj->setmTNphi(mTNphi, mNTheta+1);
-  delete [] mTNphi;
-
-  double* mTheta= new double[mNTheta+1];
-  is.GetR8s(mTheta, mNTheta+1);
-  dobj->setmTheta(mTheta, mNTheta+1);
-  delete [] mTheta;
-
-  T* mPixels= new T[mNPix];
-  PIOSReadArray(is, mPixels, mNPix);
-  dobj->setDataBlock(mPixels, mNPix);
-  delete [] mPixels;
-}
-
-template <class T>
-void FIO_SphereThetaPhi<T>::WriteSelf(POutPersist& os) const
-{
-  cout << " FIO_SphereThetaPhi:: WriteSelf " << endl;
-
-  if(dobj == NULL) 
-    {
-      cout << " WriteSelf:: dobj= null " << endl;
-      return;
-    }
-
-  char strg[256];
-  int mNTheta= dobj->SizeIndex();
-  int mNPix  = dobj->NbPixels();
- 
-  if(dobj->ptrInfo()) 
-    {
-      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d HasInfo",mNTheta,mNPix);
-      os.PutLine(strg);
-      os << dobj->Info();
-    }
-  else 
-    { 
-      sprintf(strg,"SphereThetaPhi: NSlices=%6d  NPix=%9d ",mNTheta,mNPix);
-      os.PutLine(strg);  
-    }
-
-  os.PutI4(mNTheta);  
-  os.PutI4(mNPix);
-  os.PutR8(dobj->PixSolAngle(0));
-  os.PutI4s(dobj->getmNPhi() , mNTheta);
-  os.PutI4s(dobj->getmTNphi(), mNTheta+1);
-  os.PutR8s(dobj->getmTheta(), mNTheta+1);
-  //os.Put((dobj->getDataBlock())->Data(), mNPix);
-  PIOSWriteArray(os,(dobj->getDataBlock())->Data(), mNPix);
-}
