source: Sophya/trunk/ArchTOIPipe/ProcWSophya/PSB2ring.cc@ 2755

// 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;

  int i;
  // on branche les TOIs des 4 bolos (l'un peut etre "vide")
  for(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(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;
    }
}