Index: trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.cc	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.cc	(revision 1984)
@@ -0,0 +1,311 @@
+// This may look like C code, but it's really -*- C++ -*-
+
+/////////////////////////////
+//      Bolos2ring.cc       //
+/////////////////////////////
+
+/// This is a processor for ArchTOIPipe using the Sophya Library.
+/// An example of the use of it is in the file gph425_5.cc
+/// 
+/// The goal of this processor is to give Q and U rings
+// computed with either 3 timelines coming from the
+// two channels of 1 PSB and one of the channels of a second PSB,
+// or with the differences of the two channels of each PSB.
+// The first of these two cases can be used if one channel is out
+// or order, the second one is supposed to be better when all channels
+// work fine, since it fully uses the advantage of PSB's.
+//  Both methods should give the same results when the 4 channels work ok.
+
+// Authors CR/NP, rcecile@in2p3.fr, ponthieu@isn.in2p3.fr
+//////////////////////////////////////////////////////////////////////////////////////
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <vector>
+//#include <strstream>
+#include "spherehealpix.h"
+#include "tmatrix.h"
+#include "tvector.h"
+#include "vector3d.h"
+#include "fitsfile.h"
+#include "ctimer.h"
+#include "intflapack.h"
+#include "smathconst.h"
+#include "ring33.h"
+#include "Bolos2ring.h"
+
+#define UNSEEN_HEALPIX (-1.6375e30)
+ 
+// Constructor
+Bolos2ring::Bolos2ring(SphereHEALPix<r_8>* ringQ, 
+		       SphereHEALPix<r_8>* ringU,
+		       SphereHEALPix<r_8>* ringQW,
+		       SphereHEALPix<r_8>* ringUW,
+		       const vector<r_8>& table_angle,
+		       int_4 *bolos_ok,
+		       int_4 wsz)
+  : ringq(ringQ), ringu(ringU), ringqw(ringQW), ringuw(ringUW), TableFP_(table_angle), Bolos_OK(bolos_ok)
+{
+  SetWSize(wsz);
+  totsncount_ = 0;
+  
+  if( ringq->NbPixels()<1) {
+    cout << "Bolos2ring::Bolos2ring  : bad number of pixel in ringQ "
+	 << ringq->NbPixels() << endl;
+    throw ParmError("Bolos2ring::Bolos2ring  : bad number of pixel in ringQ");
+  }
+  if( ringu->NbPixels()<1) {
+    cout << "Bolos2ring::Bolos2ring  : bad number of pixel in ringU "
+	 << ringu->NbPixels()
+	 << endl;
+    throw ParmError("Bolos2ring::Bolos2ring  : bad number of pixel in ringU");
+  }
+
+  if( ringqw->NbPixels()<1) {
+    cout << "Bolos2ring::Bolos2ring  : bad number of pixel in ringQW "
+	 << ringqw->NbPixels() << endl;
+    throw ParmError("Bolos2ring::Bolos2ring  : bad number of pixel in ringQW");
+   }
+  if( ringuw->NbPixels()<1) {
+    cout << "Bolos2ring::Bolos2ring  : bad number of pixel in ringUW "
+	 << ringuw->NbPixels() << endl;
+    throw ParmError("Bolos2ring::Bolos2ring  : bad number of pixel in ringUW");
+   }
+
+  if( ringq->NbPixels() != ringu->NbPixels())
+    throw(ParmError("Bolos2ring::Bolos2ring : rings don't have the same size!!"));
+
+  Npix_ = ringq->NbPixels();
+  int nlat = ringq->SizeIndex();
+  if(nlat != ringqw->SizeIndex()) ringqw->Resize(nlat);
+  if(nlat != ringuw->SizeIndex()) ringuw->Resize(nlat);
+  cout << "RINGS of " <<ringu->NbPixels() << " pixels" << endl;
+  
+  ringq->SetPixels(0.);
+  ringu->SetPixels(0.);
+  ringqw->SetPixels(0);
+  ringuw->SetPixels(0);
+}
+
+// Destructor
+Bolos2ring::~Bolos2ring()
+{
+}
+
+void Bolos2ring::PrintStatus(ostream& os)
+{
+  os << "____________________________________________" << endl
+     << "  Bolos2ring::PrintStatus() - wsize= "<<wsize_<< endl;
+  TOIProcessor::PrintStatus(os);
+  os << "  ProcessedSampleCount=" << ProcessedSampleCount() << endl;
+  os << "____________________________________________" << endl;
+}
+ 
+void Bolos2ring::init()
+{ 
+  cout << "Bolos2ring::init" << endl;
+
+  // on branche les TOIs des 4 bolos (l'un peut etre "vide")
+  for(int i=0; i < 4; i++)
+    {
+      char str[80];
+      sprintf(str,"Bolo%d",i);
+      cout << str << endl;
+      declareInput(str);
+    }
+
+  // on branche le pointage des 2 PSB
+  for(int i=0; i < 2; i++)
+    {
+      char str[80];
+      sprintf(str,"theta_pointing%d",i);
+      cout << str << endl;
+      declareInput(str);
+      sprintf(str,"phi_pointing%d",i);
+      cout << str << endl;
+      declareInput(str);
+    }
+  
+  name = "Bolos2ring";
+}
+
+void Bolos2ring::run()
+{
+  int snb = getMinIn();
+  int sne = getMaxIn();
+
+  if (!checkInputTOIIndex(0) && !checkInputTOIIndex(1)) {
+    cerr << " Bolos2ring::run() - Input TOIs (in1, in2 and in3, in4) not connected! "
+	 << endl;
+    throw ParmError("Bolos2ring::run() Input TOIs (in1, in2 and in3, in4) not connected!");
+  }
+
+  cout << " Bolo2ring::run() SNRange=" << snb << " - " << sne << endl;
+
+  
+  for(int i=0; i<4; i++) 
+    if( !checkInputTOIIndex(i) )
+    {
+      cerr << " Bolos2ring::run() - Input TOIs not connected! ("
+	   << i << ')' << endl;
+      throw ParmError("Bolos2ring::run() - Input TOIs not connected!");
+    }
+  
+  cout << "Bolos2ring::run() - SNRange="<< snb << '-' << sne << endl;
+  
+  try {
+    
+    cout << " *** combine bolometers without computing the difference *** " << endl;    
+
+    double v0, v1, v2, v3; // values of bolometers outputs
+    double l0, b0, l1, b1; // galactic coordinates
+    uint_8 vfg0, vfg1, vfg2, vfg3;
+    double theta0, phi0, theta1, phi1;
+    int_4 pixel0, pixel1; // pour chaque PSB
+       
+    //  Matrix* sts = new Matrix[Npix_]; // matrice de calcul des param Stokes
+    //  Vector* stm = new Vector[Npix_]; // St * vecteur de mesures
+    
+    Matrix stm(Npix_, 3);
+    RingOfMatrix3x3 ring(Npix_);
+    cout << "****************************" << endl;
+    cout << "   Ring declaration done" << endl;
+    cout << "****************************" << endl;
+    
+
+    for(int i =0; i<4; i++) TableFP_[i] *= Pi/180. ;
+
+    for(int s=snb; s<= sne; s++) {
+      
+      getData(0, s, v0, vfg0);
+      getData(1, s, v1, vfg1);
+      getData(2, s, v2, vfg2);
+      getData(3, s, v3, vfg3);
+      l0 = getData(4, s);
+      b0 = getData(5, s);
+      l1 = getData(6, s);
+      b1 = getData(7, s);
+      
+      theta0 = (90. - b0)/180. * Pi;
+      phi0 = l0/180. * Pi;
+      theta1 = (90. - b1)/180. * Pi;
+      phi1 = l1/180. *Pi;
+      
+      pixel0 = ringq->PixIndexSph(theta0, phi0); // Index of the hit pixel
+      pixel1 = ringu->PixIndexSph(theta1, phi1);
+      
+      if ((Bolos_OK[0] == 0) || (vfg0 > 0)) v0 = 0.;
+      if ((Bolos_OK[1] == 0) || (vfg1 > 0)) v1 = 0.;
+      if ((Bolos_OK[2] == 0) || (vfg2 > 0)) v2 = 0.;
+      if ((Bolos_OK[3] == 0) || (vfg3 > 0)) v3 = 0.;
+       
+      // Filling pixel0
+      ring.GetElement(pixel0)(0, 0) += 0.25 * 2;  
+      ring.GetElement(pixel0)(0, 1) += 0.25 *( cos(2*TableFP_[0]) + cos(2*TableFP_[1])); 
+      ring.GetElement(pixel0)(0, 2) += 0.25 *( sin(2*TableFP_[0]) + sin(2*TableFP_[1])); 
+      
+      ring.GetElement(pixel0)(1, 0) += 0.25 *( cos(2*TableFP_[0]) + cos(2*TableFP_[1]));
+      ring.GetElement(pixel0)(1, 1) += 0.125 *( 2 + cos(4*TableFP_[0]) + cos(4*TableFP_[1]));
+      ring.GetElement(pixel0)(1, 2) += 0.125*( sin(4*TableFP_[0]) + sin(4*TableFP_[1]));
+      
+      ring.GetElement(pixel0)(2, 0) += 0.25 *( sin(2*TableFP_[0]) + sin(2*TableFP_[1]));
+      ring.GetElement(pixel0)(2, 1) += 0.125*( sin(4*TableFP_[0]) + sin(4*TableFP_[1]));
+      ring.GetElement(pixel0)(2, 2) += 0.125*( 2 - cos(4*TableFP_[0]) - cos(4*TableFP_[1]));
+     
+      stm(pixel0, 0) += 0.5*(v0 + v1);
+      stm(pixel0, 1) += 0.5*( cos(2*TableFP_[0])*v0 + cos(2*TableFP_[1])*v1);
+      stm(pixel0, 2) += 0.5*( sin(2*TableFP_[0])*v0 + sin(2*TableFP_[1])*v1);
+           
+      if((Bolos_OK[0] == 0) || (Bolos_OK[1] == 0) || (vfg0 >0) || (vfg1 >0))
+      {
+	ringqw->PixVal(pixel0) += 1 ; // nombre de hits dans le pixel
+	ringuw->PixVal(pixel0) += 1 ; 
+      }
+      else 
+      {
+	ringqw->PixVal(pixel0) += 2 ; 
+	ringuw->PixVal(pixel0) += 2 ; 
+      }
+      
+      ring.GetElement(pixel1)(0, 0) += 0.25 * 2; 
+      ring.GetElement(pixel1)(0, 1) += 0.25 *( cos(2*TableFP_[2]) + cos(2*TableFP_[3])); 
+      ring.GetElement(pixel1)(0, 2) += 0.25 *( sin(2*TableFP_[2]) + sin(2*TableFP_[3])); 
+      
+      ring.GetElement(pixel1)(1, 0) += 0.25 *( cos(2*TableFP_[2]) + cos(2*TableFP_[3]));
+      ring.GetElement(pixel1)(1, 1) += 0.125*( 2 + cos(4*TableFP_[2]) + cos(4*TableFP_[3]));
+      ring.GetElement(pixel1)(1, 2) += 0.125*( sin(4*TableFP_[2]) + sin(4*TableFP_[3])); 
+      
+      ring.GetElement(pixel1)(2, 0) += 0.25 *( sin(2*TableFP_[2]) + sin(2*TableFP_[3]));
+      ring.GetElement(pixel1)(2, 1) += 0.125*( sin(4*TableFP_[2]) + sin(4*TableFP_[3]));
+      ring.GetElement(pixel1)(2, 2) += 0.125*( 2 - cos(4*TableFP_[2]) - cos(4*TableFP_[3]));
+      
+      stm(pixel1, 0) += 0.5*( v2 + v3);
+      stm(pixel1, 1) += 0.5*( cos(2*TableFP_[2])*v2 + cos(2*TableFP_[3])*v3);
+      stm(pixel1, 2) += 0.5*( sin(2*TableFP_[2])*v2 + sin(2*TableFP_[3])*v3);
+      
+      if((Bolos_OK[2] == 0) || (Bolos_OK[3] == 0)|| (vfg2 >0) || (vfg3 >0)) 
+      {  	
+	ringqw->PixVal(pixel1) += 1 ; // nombre de hits dans le pixel
+	ringuw->PixVal(pixel1) += 1 ;
+      }
+      else 
+      {
+	ringqw->PixVal(pixel1) += 2 ; 
+	ringuw->PixVal(pixel1) += 2 ; 
+      }
+      
+      
+      if (!(s%10000)) cout << " sample " << s << " ****** filling matrices  ****** " << endl;
+      
+    }// fin de la boucle sur les samplenums qui normalement a rempli
+    // la matrice et le vecteur projete. Y'a plus qu'a inverser...
+    
+    // Ca devrait suffire de definir un ringw seulement dans ce cas la,
+    // mais je garde le ringqw pour pas bugguer l'entree.
+    
+    cout << endl << "Start to fill rings :" << endl;
+	   
+    for(int pix=0 ; pix < Npix_ ; pix++) {// pour chaque pixel
+
+      if( ringqw->PixVal(pix) < 4.0 ) 
+	//pixel pas vu, 4 a cause des 2 PSB ici, a changer
+	ringq->PixVal(pix) = ringu->PixVal(pix) = UNSEEN_HEALPIX; 
+      
+      else{ 
+      
+	Matrix3x3 sts;
+	Matrix stsm1(3, 3);
+	Vector s(3);
+	Vector stokes_param(3);
+	
+	s(0) = stm(pix, 0);
+	s(1) = stm(pix, 1);
+	s(2) = stm(pix, 2);
+	sts = ring.GetElement(pix); // get the sts matrix of pix
+	sts = sts.Inverse();
+
+	// voir si on peut arreter de bricoler pour que Matrix3x3 soit aussi une Matrix
+	for(int r=0; r<=2; r++)
+	  for(int c=0; c<=2; c++)
+	    stsm1(r, c) = sts(r, c);
+
+	stokes_param = stsm1 * s;	
+  	ringq->PixVal(pix) = stokes_param(0);
+  	ringu->PixVal(pix) = stokes_param(2);
+	
+      }
+      
+    }
+    
+  }
+  
+  
+
+  catch( PThrowable& exc )
+  {
+    cerr << " Exception: " << exc.Msg() << endl;
+  }
+  
+  return;
+} 
Index: trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.h
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.h	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/Bolos2ring.h	(revision 1984)
@@ -0,0 +1,61 @@
+// This may look like C code, but it's really -*- C++ -*-
+
+#if !defined(_Bolos2RING_H_SEEN_)
+#define _Bolos2RING_H_SEEN_
+
+#include <iostream.h>
+#include <string>
+#include "toiprocessor.h"
+#include "spherehealpix.h"
+#include "localmap.h"
+#include "flagtoidef.h"
+
+class Bolos2ring : public TOIProcessor
+{
+public:
+  Bolos2ring(SphereHEALPix<r_8>* ringQ, 
+	     SphereHEALPix<r_8>* ringU,
+	     SphereHEALPix<r_8>* ringQW,
+	     SphereHEALPix<r_8>* ringUW,
+	     const vector<r_8>& table_angle,
+	     int_4 *Bolos_OK,
+	     int_4 wsz=8192);
+
+  virtual ~Bolos2ring();
+  inline void SetWSize(int_4 wsz) { wsize_ = (wsz<4) ? 4 : wsz; }
+  inline int_8 ProcessedSampleCount() const { return totsncount_; }
+  void SetNSide(int_4 nside);
+  virtual void init();  virtual void run();
+  
+  virtual void PrintStatus(ostream& os);
+
+protected:
+  bool Bolos2ringdone;
+
+  SphereHEALPix<r_8>* ringq;   // carte de Q
+  SphereHEALPix<r_8>* ringu;   // carte de U
+  
+
+  SphereHEALPix<r_8>* ringqw;   // carte du nb de hit
+  SphereHEALPix<r_8>* ringuw;
+  vector<r_8> TableFP_;
+  int_4* Bolos_OK;
+
+  int wsize_;
+  int_8 totsncount_;
+  int_8 Npix_; 
+  
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/ArchTOIPipe/ProcWSophya/PSB2ring.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/PSB2ring.cc	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/PSB2ring.cc	(revision 1984)
@@ -0,0 +1,247 @@
+// This may look like C code, but it's really -*- C++ -*-
+
+/////////////////////////////
+//      PSB2ring.cc       //
+/////////////////////////////
+
+/// This is a processor for ArchTOIPipe using the Sophya Library.
+/// An example of the use of it is in the file gph425_5.cc
+/// 
+/// The goal of this processor is to give Q and U rings
+// computed with either 3 timelines coming from the
+// two channels of 1 PSB and one of the channels of a second PSB,
+// or with the differences of the two channels of each PSB.
+// The first of these two cases can be used if one channel is out
+// or order, the second one is supposed to be better when all channels
+// work fine, since it fully uses the advantage of PSB's.
+//  Both methods should give the same results when the 4 channels work ok.
+
+// Authors CR/NP, rcecile@in2p3.fr, ponthieu@isn.in2p3.fr
+//////////////////////////////////////////////////////////////////////////////////////
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <vector>
+//#include <strstream>
+#include "tmatrix.h"
+#include "tvector.h"
+#include "vector3d.h"
+#include "ctimer.h"
+#include "PSB2ring.h"
+
+#define UNSEEN_HEALPIX (-1.6375e30)
+
+
+// Constructor
+PSB2ring::PSB2ring(SphereHEALPix<r_8>* ringQ, 
+		   SphereHEALPix<r_8>* ringU,
+		   SphereHEALPix<r_8>* ringQW,
+		   SphereHEALPix<r_8>* ringUW,
+		   const vector<r_8>& table_angle,
+		   int_4 wsz)
+  : ringq(ringQ), ringu(ringU), ringqw(ringQW), ringuw(ringUW), TableFP_(table_angle)
+{
+  SetWSize(wsz);
+  totsncount_ = 0;
+  
+  if( ringq->NbPixels()<1) {
+    cout << "PSB2ring::PSB2ring  : bad number of pixel in ringQ "
+	 << ringq->NbPixels() << endl;
+    throw ParmError("PSB2ring::PSB2ring  : bad number of pixel in ringQ");
+  }
+  if( ringu->NbPixels()<1) {
+    cout << "PSB2ring::PSB2ring  : bad number of pixel in ringU "
+	 << ringu->NbPixels()
+	 << endl;
+    throw ParmError("PSB2ring::PSB2ring  : bad number of pixel in ringU");
+  }
+
+  if( ringqw->NbPixels()<1) {
+    cout << "PSB2ring::PSB2ring  : bad number of pixel in ringQW "
+	 << ringqw->NbPixels() << endl;
+    throw ParmError("PSB2ring::PSB2ring  : bad number of pixel in ringQW");
+   }
+  if( ringuw->NbPixels()<1) {
+    cout << "PSB2ring::PSB2ring  : bad number of pixel in ringUW "
+	 << ringuw->NbPixels() << endl;
+    throw ParmError("PSB2ring::PSB2ring  : bad number of pixel in ringUW");
+   }
+
+  if( ringq->NbPixels() != ringu->NbPixels())
+    throw(ParmError("PSB2ring::PSB2ring : rings don't have the same size!!"));
+
+  Npix_ = ringq->NbPixels();
+  int nlat = ringq->SizeIndex();
+  if(nlat != ringqw->SizeIndex()) ringqw->Resize(nlat);
+  if(nlat != ringuw->SizeIndex()) ringuw->Resize(nlat);
+  cout << "RINGS of " <<ringu->NbPixels() << " pixels" << endl;
+  
+  ringq->SetPixels(0.);
+  ringu->SetPixels(0.);
+  ringqw->SetPixels(0);
+  ringuw->SetPixels(0);
+}
+
+// Destructor
+PSB2ring::~PSB2ring()
+{
+}
+
+void PSB2ring::PrintStatus(ostream& os)
+{
+  os << "____________________________________________" << endl
+     << "  PSB2ring::PrintStatus() - wsize= "<<wsize_<< endl;
+  TOIProcessor::PrintStatus(os);
+  os << "  ProcessedSampleCount=" << ProcessedSampleCount() << endl;
+  os << "____________________________________________" << endl;
+}
+ 
+void PSB2ring::init()
+{ 
+  cout << "PSB2ring::init" << endl;
+
+  // on branche les TOIs des 4 bolos (l'un peut etre "vide")
+  for(int i=0; i < 4; i++)// if (PSB_ok[i] > 0)
+    {
+      char str[80];
+      sprintf(str,"Bolo%d",i);
+      cout << str << endl;
+      declareInput(str);
+    }
+
+  // on branche le poinatge des 2 PSB
+  for(int i=0; i < 2; i++)// if (PSB_ok[i] > 0)
+    {
+      char str[80];
+      sprintf(str,"theta_pointing%d",i);
+      cout << str << endl;
+      declareInput(str);
+      sprintf(str,"phi_pointing%d",i);
+      cout << str << endl;
+      declareInput(str);
+    }
+
+  declareOutput("out_toiQ");
+  declareOutput("out_toiU");
+ 
+  name = "PSB2ring";
+}
+
+void PSB2ring::run()
+{
+  int snb = getMinIn();
+  int sne = getMaxIn();
+
+  if (!checkInputTOIIndex(0) && !checkInputTOIIndex(1)) {
+    cerr << " PSB2ring::run() - Input TOIs (in1, in2 and in3, in4) not connected! "
+	 << endl;
+    throw ParmError("PSB2ring::run() Input TOIs (in1, in2 and in3, in4) not connected!");
+  }
+
+  cout << " Bolo2ring::run() SNRange=" << snb << " - " << sne << endl;
+  
+  for(int i=0; i<4; i++) //if (PSB_ok[i/2,i%2])
+    if( !checkInputTOIIndex(i) )
+    {
+      cerr << " PSB2ring::run() - Input TOIs not connected! ("
+	   << i << ')' << endl;
+      throw ParmError("PSB2ring::run() - Input TOIs not connected!");
+    }
+  
+  // pas de sortie a verifier
+  
+  cout << "PSB2ring::run() - SNRange="<< snb << '-' << sne << endl;
+  
+  try {
+
+    Timer tm("PSB2ring::run()");
+    cout << "****************************" << endl;
+    cout << "****************************" << endl;  
+    cout << "4 bolos working and ordered" << endl;
+    
+    // Generation des timelines
+    double v0,v1,v2,v3; // values of bolometers outputs
+    double l0, b0, l1, b1; // galactic coordinates
+    uint_8 vfg0,vfg1,vfg2,vfg3;
+    double Q,U;    
+    
+    for(int s=snb;s<=sne;s++) {
+
+      getData(0, s, v0, vfg0);
+      getData(1, s, v1, vfg1);
+      getData(2, s, v2, vfg2);
+      getData(3, s, v3, vfg3);
+
+      l0 = getData(4, s);
+      b0 = getData(5, s);
+      l1 = getData(6, s);
+      b1 = getData(7, s);
+
+      double theta0 = (90. - b0)/180 * 3.14159;
+      double phi0 = l0/180. * 3.14159;
+      double theta1 = (90. - b1)/180 * 3.14159;
+      double phi1 = l1/180 *3.14159;
+
+
+
+      //     cout << " ********** Q(t) and U(t) finished *********** " << endl;
+ 
+      //    cout << " ********** Projection on the ring ********* " << endl;
+      
+      //Projection
+      
+      if((vfg0 == 0) && (vfg1 == 0))
+      {
+	Q = v0 - v1;
+	int_4 pixelQ = ringq->PixIndexSph(theta0, phi0); // numero du pixel touche
+	ringqw->PixVal(pixelQ) += 1 ; // nombre de hits dans le pixel
+	ringq->PixVal(pixelQ)  += Q ;
+      }
+      
+      if((vfg2 == 0) && (vfg3 == 0))
+      {
+	U = v2 - v3;      
+	int_4 pixelU = ringu->PixIndexSph(theta1, phi1);
+	ringuw->PixVal(pixelU) += 1 ;      
+	ringu->PixVal(pixelU)  += U ;
+      }
+      
+      if (!(s%10000)) cout << s << " ****** filling rings  ****** " << endl;
+	
+      putData(0, s, Q);   
+      putData(1, s, U); 
+    }
+
+    // Moyennage des pixels
+    // comme on bosse directement avec des timelines de Q et U
+    // dans ce cas precis, on a juste a prendre la moyenne des valeurs
+    // qui tombent dans un pixel pour avoir la carte
+    // on ne doit inverser de matrice etc que dans le cas
+    // trois bolos.
+
+    for(int pix=0 ; pix < Npix_ ; pix++) {
+
+      if( ringqw->PixVal(pix) < 1.0 ) //pixel pas vu
+	ringq->PixVal(pix) = ringu->PixVal(pix) = UNSEEN_HEALPIX; 
+      else
+	ringq->PixVal(pix) = ringq->PixVal(pix) / ringqw->PixVal(pix);
+      
+      if( ringuw->PixVal(pix) < 1.0 )
+	ringu->PixVal(pix) = ringu->PixVal(pix) = UNSEEN_HEALPIX; 
+      else
+	ringu->PixVal(pix) = ringu->PixVal(pix) / ringuw->PixVal(pix);      
+
+    }
+    cout << " ********** rings filled  ********* " << endl;
+  }
+  
+
+  catch( PThrowable& exc )
+    {
+      cerr << " Exception: " << exc.Msg() << endl;
+    }
+}
+
+
+
Index: trunk/ArchTOIPipe/ProcWSophya/PSB2ring.h
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/PSB2ring.h	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/PSB2ring.h	(revision 1984)
@@ -0,0 +1,64 @@
+// This may look like C code, but it's really -*- C++ -*-
+
+#if !defined(_PSB2RING_H_SEEN_)
+#define _PSB2RING_H_SEEN_
+
+#include <iostream.h>
+#include <string>
+#include "toiprocessor.h"
+#include "spherehealpix.h"
+#include "localmap.h"
+#include "flagtoidef.h"
+
+class PSB2ring : public TOIProcessor
+{
+public:
+  PSB2ring(SphereHEALPix<r_8>* ringQ, 
+	   SphereHEALPix<r_8>* ringU,
+	   SphereHEALPix<r_8>* ringQW,
+	   SphereHEALPix<r_8>* ringUW,
+	   const vector<r_8>& table_angle,
+	   int_4 wsz=819);
+
+  virtual ~PSB2ring();
+  inline void SetWSize(int_4 wsz) { wsize_ = (wsz<4) ? 4 : wsz; }
+  inline int_8 ProcessedSampleCount() const { return totsncount_; }
+  void SetNSide(int_4 nside);
+  virtual void init();
+  virtual void run();
+  
+  virtual void PrintStatus(ostream& os);
+
+protected:
+  bool PSB2ringdone;
+
+  SphereHEALPix<r_8>* ringq;   // carte de Q
+  SphereHEALPix<r_8>* ringu;   // carte de U
+  
+  string ringqfile;
+  string ringufile;
+
+  vector<r_8> TableFP_;
+  SphereHEALPix<r_8>* ringqw;   // carte du nb de hit
+  SphereHEALPix<r_8>* ringuw;
+  
+  string ringqwfile;
+  string ringuwfile;
+
+  int wsize_;
+  int_8 totsncount_;
+  int_8 Npix_;
+};
+
+#endif
+
+
+
+
+
+
+
+
+
+
+
Index: trunk/ArchTOIPipe/ProcWSophya/ring33.cc
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/ring33.cc	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/ring33.cc	(revision 1984)
@@ -0,0 +1,22 @@
+#include "machdefs.h"
+#include "pexceptions.h"
+#include "ring33.h"
+
+
+// -------------------------------------------------------
+// Classe simple de matrice 3x3 - sans allocation memoire
+// -------------------------------------------------------
+
+// constructor
+RingOfMatrix3x3::RingOfMatrix3x3(int_4 npix)
+{
+  if (npix < 1) throw RangeCheckError("RingOfMatrix3x3::RingMatrix3x3() size < 1 !");
+  size_ = npix;
+  mtxtab_ = new Matrix3x3[npix];
+}
+
+// Destructor
+RingOfMatrix3x3::~RingOfMatrix3x3()
+{
+  delete[] mtxtab_;
+}
Index: trunk/ArchTOIPipe/ProcWSophya/ring33.h
===================================================================
--- trunk/ArchTOIPipe/ProcWSophya/ring33.h	(revision 1984)
+++ trunk/ArchTOIPipe/ProcWSophya/ring33.h	(revision 1984)
@@ -0,0 +1,147 @@
+#include "machdefs.h"
+#include "pexceptions.h"
+
+/*
+#include <math.h>
+#include <iostream.h>
+#include <fstream.h>
+
+#include "tarrinit.h"
+#include "array.h"
+#include "ctimer.h"
+#include "timing.h"
+*/
+
+
+
+// -------------------------------------------------------
+// Classe simple de matrice 3x3 - sans allocation memoire
+// -------------------------------------------------------
+
+class Matrix3x3 {
+public:
+  inline        Matrix3x3(bool raz=true) 
+  { 
+    if (raz) 
+      for(int i=0; i<9; i++) a_[i]=0.; 
+  } 
+  inline	Matrix3x3(const Matrix3x3 & m)
+  { 
+    for(int i=0; i<9; i++)  {a_[i] = m.a_[i]; am1_[i] = m.am1_[i];}
+  }
+
+  inline        ~Matrix3x3() { }
+
+  inline Matrix3x3& Set(const Matrix3x3 & m)
+                   { for(int i=0; i<9; i++) a_[i] = m.a_[i]; return (*this); }
+  inline Matrix3x3& operator = (const Matrix3x3 & a)
+                   { return Set(a); }
+
+  inline double operator()(int r, int c) const 
+                { return a_[r*3+c]; }
+  inline double& operator()(int r, int c) 
+                { return a_[r*3+c]; }
+
+  inline int NRows() const {return 3; }
+  inline int NCols() const {return 3; }
+
+  inline Matrix3x3& AddElt(const Matrix3x3 & m)
+                   { for(int i=0; i<9; i++) a_[i] += m.a_[i]; return (*this); }
+  inline Matrix3x3& MulElt(const Matrix3x3 & m)
+                   { for(int i=0; i<9; i++) a_[i] *= m.a_[i]; return (*this); }
+  
+  // Multiplication de deux Matrix3x3 entre elles  
+/*    inline Matrix3x3& MulMat(const Matrix3x3 & m) */
+/*      { */
+/*        Matrix3x3 c; */
+/*        // c.Set(0.); */
+/*        for(int r=0; r<=2; r++) */
+/*  	for(int c=0; c<=2; c++) */
+/*  	  c.a_[r*3+c] = 0.; */
+      
+/*        for(int r=0; r<=2; r++) */
+/*  	for(int c=0; c<=2; c++) */
+/*  	  for(int k=0; k<=2; k++) */
+/*  	    c.a_[3*r+c] += a_[3*r+k] * m.a_[3*k+c]; */
+      
+/*        a_ = c; */
+/*        return(*this); */
+/*      } */
+
+  inline Matrix3x3& Transpose()
+    {
+      double c;
+      c = a_[1]; a_[1] = a_[3]; a_[3] = c;
+      c = a_[2]; a_[2] = a_[6]; a_[6] = c;
+      c = a_[5]; a_[5] = a_[7]; a_[7] = c;
+      return (*this);
+    }
+
+  inline Matrix3x3& Inverse()
+    {
+      // Calcul du determinant global
+      double determ;
+      determ =  a_[0]*a_[4]*a_[8] + a_[3]*a_[7]*a_[2] + a_[6]*a_[1]*a_[5];
+      determ -= a_[2]*a_[4]*a_[6] + a_[5]*a_[7]*a_[0] + a_[8]*a_[1]*a_[3];
+      // prevoir un break si determ <= threshold
+
+      // Calcul des cofacteurs
+      am1_[0] =  (a_[4]*a_[8] - a_[5]*a_[7])/determ ; 
+      am1_[1] = -(a_[3]*a_[8] - a_[5]*a_[6])/determ ; 
+      am1_[2] =  (a_[3]*a_[7] - a_[4]*a_[6])/determ ; 
+      am1_[3] = -(a_[1]*a_[8] - a_[7]*a_[2])/determ ; 
+      am1_[4] =  (a_[0]*a_[8] - a_[6]*a_[2])/determ ; 
+      am1_[5] = -(a_[0]*a_[7] - a_[6]*a_[1])/determ ; 
+      am1_[6] =  (a_[1]*a_[5] - a_[4]*a_[2])/determ ; 
+      am1_[7] = -(a_[0]*a_[5] - a_[3]*a_[2])/determ ; 
+      am1_[8] =  (a_[0]*a_[4] - a_[3]*a_[1])/determ ; 
+      /****************************************/
+      /* c'est la version qui marche, mais la convention semble 
+         differente de celle de [3*r + c]...? */
+
+
+/*        // Transposition */
+/*        am1_ = am1_.Transpose() ; */
+/*        a_ = am1_; // marche pas bien : "request for member `Transpose' in `Matrix3x3::am1_', */
+/*                   //which is of non-aggregate type 'double[9]' */
+
+      // Calcul de l'inverse avec transpostion incluse
+/*        am1_[0] =  (a_[4]*a_[8] - a_[5]*a_[7])/determ ; */
+/*        am1_[3] = -(a_[3]*a_[8] - a_[5]*a_[6])/determ ; */
+/*        am1_[6] =  (a_[3]*a_[7] - a_[4]*a_[6])/determ ; */
+/*        am1_[1] = -(a_[1]*a_[8] - a_[7]*a_[2])/determ ; */
+/*        am1_[4] =  (a_[0]*a_[8] - a_[6]*a_[2])/determ ; */
+/*        am1_[7] = -(a_[0]*a_[7] - a_[6]*a_[1])/determ ; */
+/*        am1_[2] =  (a_[1]*a_[5] - a_[4]*a_[2])/determ ; */
+/*        am1_[5] = -(a_[0]*a_[5] - a_[3]*a_[2])/determ ; */
+/*        am1_[8] =  (a_[0]*a_[4] - a_[3]*a_[1])/determ ; */
+      a_ = am1_;
+      return (*this);}
+
+
+protected:
+  double a_[9], am1_[9];
+  
+};
+
+
+class RingOfMatrix3x3 {
+public:
+                RingOfMatrix3x3(int_4 npix);
+  virtual      ~RingOfMatrix3x3();
+
+  inline int_4  Size() { return size_; }
+  inline Matrix3x3& GetElement(int_4 n) 
+  {
+    if ((n<0) || (n>=size_))
+      throw RangeCheckError("RingOfMatrix3x3::GetElement() - Out of range index");
+    return mtxtab_[n];
+  }
+  inline Matrix3x3& operator()(int_4 n) { return GetElement(n); }
+
+protected:
+  Matrix3x3 *mtxtab_;
+  int_4 size_;  
+};
+
+
