source: Sophya/trunk/ArchTOIPipe/ProcWSophya/toi2map.cc@ 2390

// ArchTOIPipe           (C)     CEA/DAPNIA/SPP IN2P3/LAL
//                               Eric Aubourg
//                               Christophe Magneville
//                               Reza Ansari
// $Id: toi2map.cc,v 1.33 2003-01-27 07:56:04 aubourg Exp $

#include "machdefs.h"
#include "toimanager.h"
#include "pexceptions.h"
#include "ctimer.h"
#include "toi2map.h"
// La valeur "Pi" doit etre celle de smathconst.h a cause du test sur theta
#include "smathconst.h"

#include "fitsspherehealpix.h"


////////////////////////////////////////////////////////////////////////
TOI2Map::TOI2Map(string mapFName, string wmapFName, int nside) 
  : mMap(new SphereHEALPix<r_8>(nside)), mWMap(new SphereHEALPix<r_8>(nside)),
    mWMapInternal(false), totnscount(0), usednscount(0),
    mMapFName(mapFName), mWMapFName(wmapFName)
{
  commonConst();
}

////////////////////////////////////////////////////////////////////////
TOI2Map::TOI2Map(PixelMap<r_8>* map,PixelMap<r_8>* wmap)
  : mMap(map), mWMap(wmap), mWMapInternal(false), totnscount(0), usednscount(0),
    mMapFName(""), mWMapFName("")
{
  commonConst();
}

void TOI2Map::commonConst() {
 SetEquinox();
 SetCoorIn();
 SetCoorMap();
 SetCalibrationFactor();
 SetTestFlag();
 SetTestMin();
 SetTestMax();

 if(mMap->NbPixels()<1) {
  cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mMap "
      <<mMap->NbPixels()<<endl;
  throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
 }
 mMap->SetPixels(0.);
 int nlat = mMap->SizeIndex();
 string typmap = mMap->TypeOfMap();

 if(mWMap==NULL) {
   // We would need a cloning function in maps. $$TBD$$ $CHECK$
   if (typmap == "LOCAL") {
     mWMap = new LocalMap<r_8>(*(LocalMap<r_8>*)mMap,false);
   } else if (typmap == "RING") {
     mWMap = new SphereHEALPix<r_8>(nlat);
   } else {
     cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
     throw ParmError("TOI2Map::TOI2Map() - bad type of map");
   }
   mWMapInternal = true;
 } else {
   mWMapInternal = false;
   if(nlat != mWMap->SizeIndex()) {
     cout<<"TOI2Map::TOI2Map() Bad size for sphere mWMap, does not "
         <<"correspond to mMap : "<<mMap->SizeIndex()<<", "<<mWMap->SizeIndex()<<endl;
     if (typmap == "LOCAL") {
       throw ParmError("TOI2Map::TOI2Map() - Different sizes for map and wmap");
     } else if(typmap == "RING") {
       ((SphereHEALPix<r_8> *)mWMap)->Resize(nlat);
       cout<<"Resize have been done..."<<endl;
     } else {
       cout << "TOI2Map::TOI2Map() cannot handle map of type " << typmap << endl;
       throw ParmError("TOI2Map::TOI2Map() - bad type of map");
     }
   }
 }
 if(mWMap->NbPixels()<1) {
   cout<<"TOI2Map::TOI2Map() Bad number of pixels in sphere mWMap "
       <<mWMap->NbPixels()<<endl;
   throw ParmError("TOI2Map::TOI2Map() - Bad number of pixels in sphere");
 }
 mWMap->SetPixels(0);
}

TOI2Map::~TOI2Map()
{
 if(mWMap && mWMapInternal) delete mWMap;
}

////////////////////////////////////////////////////////////////////////
void TOI2Map::Print(::ostream & os)
{
 os<<"TOI2Map::Print -- Map type " << mMap->TypeOfMap() << " SizeIndex = "<<mMap->SizeIndex()<<endl
   <<"   - Equinoxe="<<mActualYear<<endl
   <<"   - TypCoorIn:  "<<mTypCoorIn<<" = "<<DecodeTypAstro(mTypCoorIn)<<endl
   <<"   - TypCoorMap: "<<mTypCoorMap<<" = "<<DecodeTypAstro(mTypCoorMap)<<endl
   <<"  - Tests: Flag("<<mTFlag<<") bad="<<mBadFlag
   <<" / Value Min("<<mTMin<<")="<<mValMin
   <<" , Max("<<mTMax<<")="<<mValMax<<endl;
 os<<" TotalProcessedSamples="<<ProcessedSampleCount()
   <<" UsedSampleCount="<<UsedSampleCount()<<endl;
}

////////////////////////////////////////////////////////////////////////
void TOI2Map::init() {
  cout << "TOI2Map::init" << endl;
  declareInput("Coord1In");     // input index 0
  declareInput("Coord2In");     // input index 1
  declareInput("BoloIn");       // input index 2
}

////////////////////////////////////////////////////////////////////////
// define SANS_BUFFER
void TOI2Map::run()
{
long snb = getMinIn();
long sne = getMaxIn();

if(snb>sne) {
  cout<<"TOI2Map::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
  throw ParmError("TOI2Map::run() - Bad sample interval");
}
if(!checkInputTOIIndex(0) || !checkInputTOIIndex(1) || !checkInputTOIIndex(2)) {
  cout<<"TOI2Map::run() - Input TOI (Coord1In or Coord2In or BoloIn) not connected! "<<endl;
  throw ParmError("TOI2Map::run() Output TOI (Coord1In or Coord2In or BoloIn) not connected!");
}
if( !(mTypCoorIn&TypCoordEq || mTypCoorIn&TypCoordGal) ) {
  cout<<"TOI2Map::run() - Input Coordinates not Eq or Gal! "<<endl;
  throw ParmError("TOI2Map::run() - Input Coordinates not Eq or Gal!");
}
if( !(mTypCoorMap&TypCoordEq || mTypCoorMap&TypCoordGal) ) {
  cout<<"TOI2Map::run() - Output Coordinates not Eq or Gal! "<<endl;
  throw ParmError("TOI2Map::run() - Output Coordinates not Eq or Gal!");
}

//---------------------------------------------------------
#define NFILL 25
try {

int ii;
uint_4 mNSnFill=0, mNpixFill=0, NFill[NFILL], BadCoorRange=0, nbadsn=0;
for(ii=0;ii<NFILL;ii++) NFill[ii]=0;
double mjd = MJDfrYear(mActualYear);

cout<<"TOI2Map::run() from "<<snb<<" to "<<sne;
#ifndef SANS_BUFFER
int bufsz = 100;
uint_8* bfgbolo = new uint_8[bufsz];
double* bbolo   = new double[bufsz];
double* bc1     = new double[bufsz];
double* bc2     = new double[bufsz];
int i0 = -1;
cout<<" (getData() bufferized)"<<endl;
#else
cout<<" (getData() not bufferized)"<<endl;
#endif
 cout << "calib = " << mCalibFactor << endl;
// Remplissage des spheres
for(int s=snb;s<=sne;s++) {
  uint_8 fgbolo = 0;
  double bolo,coord1,coord2;
  //              Equatoriales   /   Galactiques
  // coord1,2 =   alpha,delta    /   gLon,gLat

#ifndef SANS_BUFFER
  if(i0<0 || s<i0 || s>=i0+bufsz) {
    i0 = s;
    int nget = (sne-s+1<bufsz)? nget=sne-s+1: bufsz;
    getData(0,i0,nget,bc1);
    getData(1,i0,nget,bc2);
    getData(2,i0,nget,bbolo,bfgbolo);
    totnscount += nget;
  }
  bolo = bbolo[s-i0]; fgbolo = bfgbolo[s-i0];
  coord1 = bc1[s-i0]; coord2 = bc2[s-i0];
#else
  getData(2,s,bolo,fgbolo);
  coord1 = getData(0,s);
  coord2 = getData(1,s);
  totnscount++;
#endif

  // Comme il n'y a pas de toi en sortie, il faut dire
  // aux processeur/toi que l'on a plus besoin des donnees.
  if (s%100 == 0) wontNeedBefore(s-1);
  //if(fgbolo>0) cout<<fgbolo<<" , "<<(fgbolo&mBadFlag)<<endl;
  if(mTFlag && (fgbolo&mBadFlag)) {nbadsn++; continue;}
  if(mTMin && (bolo<mValMin))     {nbadsn++; continue;}
  if(mTMax && (bolo>mValMax))     {nbadsn++; continue;}
  usednscount++;

  // sphere phi   entre [0,2*Pi] en radians
  // sphere theta entre [0,Pi]   en radians
  double phi=-1.;
  CoordConvertToStd(mTypCoorIn,&coord1,&coord2);

  if(mTypCoorIn&TypCoordEq && mTypCoorMap&TypCoordGal) { // Eq -> Gal
    EqtoGal(mjd,coord1,coord2,&coord1,&coord2);
    phi   = coord1 * Pi/180.;
  } else if(mTypCoorIn&TypCoordGal && mTypCoorMap&TypCoordEq) { // Gal -> Eq
    GaltoEq(mjd,coord1,coord2,&coord1,&coord2);
    phi   = coord1 * Pi/12.;
  } else if(mTypCoorMap&TypCoordGal) { // Gal -> Gal
    phi   = coord1 * Pi/180.;
  } else if(mTypCoorMap&TypCoordEq) { // Eq -> Eq
    phi   = coord1 * Pi/12.;
  }
  ToCoLat(&coord2,TypUniteD);
  double theta = coord2 * Pi/180.;

  if(phi<0. | phi>=2*Pi || theta<0. || theta>Pi)
                     {BadCoorRange++; continue;}

  int_4 ipix = mMap->PixIndexSph(theta,phi);
  if ((ipix < 0) || (ipix >= mMap->NbPixels()) ) continue;
  /*
    if (mNSnFill % 1000 == 0) {
    cout << ipix << " " << (*mMap)(ipix) << " + " << bolo*mCalibFactor <<
    " " << ((*mWMap)(ipix)) << endl;
      }
  */
  (*mMap)(ipix) += bolo*mCalibFactor;    
  ((*mWMap)(ipix)) += 1;
  mNSnFill++;
}

 cout<<"TOI2Map::run(): End of SN loop - TotalProcessedSamples="
     <<ProcessedSampleCount()<<" UsedSampleCount="<<UsedSampleCount()<<endl;

// Remplissage des spheres
 for(int_4 i=0;i<mMap->NbPixels();i++) {
   r_8 wf = (*mWMap)(i);
   /*
     if (i%100 == 0) {
     cout << "pix " << i << " w = " << wf << " map = " <<
     (*mMap)(i);
     }
   */
   if(wf>0.) {mNpixFill++; (*mMap)(i) /= wf;}
   /*
     if (i%100 == 0) {
     cout << " -> " << (*mMap)(i);
     }
   */
   int_4 nf = int_4(wf);
   if(nf>=NFILL) nf=NFILL-1; NFill[nf]++;
 }

 cout<<"TOI2Map::run(): mNpixTot="<<mMap->NbPixels()
     <<"  mNpixFill="<<mNpixFill
     <<"  --> FracSky="<<mNpixFill*100./(double)mMap->NbPixels()<<"%"
     <<"  NFill["<<NFILL<<"] ="<<endl;
 for(ii=0;ii<NFILL;ii++) {cout<<NFill[ii]<<" "; if(ii%10==9) cout<<endl;}
 cout<<endl;
 cout<<"  mNSnFill="<<mNSnFill<<" NBadSn="<<nbadsn
     <<"  BadCoorRange="<<BadCoorRange<<endl;

#ifndef SANS_BUFFER
delete [] bbolo; delete [] bfgbolo;
delete [] bc1;   delete [] bc2;
#endif

/* if (mMapFName != "") {
   cout << "saving maps" << endl;
   {
     FitsOutFile sfits(mMapFName,FitsFile::clear);
     sfits << * (SphereHEALPix<r_8>*)mMap;
   }
   {
     FitsOutFile sfits(mWMapFName,FitsFile::clear);
     sfits << * (SphereHEALPix<r_8>*)mWMap;
   }
 } 
*/
//---------------------------------------------------------
} catch (PException & exc) {
  cout<<"TOI2Map: Catched Exception "<<(string)typeid(exc).name()
      <<"\n .... Msg= "<<exc.Msg()<<endl;
}

return;                                                                            
}