source: Sophya/trunk/FrEROS/AnaLC/lcproc.cc@ 3324

#include "lcproc.h"

#include "sopnamsp.h"

LightCurveProc::LightCurveProc(string outppf, int dtblks, int dtblkti)
  : dts(dtblks) , dtti(dtblkti), _outname(outppf)
{
  // Defining column names for the dts : One row / star 
  dts.AddFloatColumn("NumEt");
  dts.AddFloatColumn("FluxRef");
  dts.AddFloatColumn("XPos");
  dts.AddFloatColumn("YPos");
  dts.AddFloatColumn("nmes");
  dts.AddFloatColumn("MeanFlux");
  dts.AddFloatColumn("SigFlux");


  // Defining column names for the dtti : One row / exposure ( = TimeInfo_ 
  dtti.AddFloatColumn("TStart");
  dtti.AddFloatColumn("Expose");
  dtti.AddFloatColumn("Fond");
  dtti.AddFloatColumn("SigFond");
  dtti.AddFloatColumn("SigX");
  dtti.AddFloatColumn("SigY");
  dtti.AddFloatColumn("Absorption");
  dtti.AddIntegerColumn("FgBad");
}

LightCurveProc::~LightCurveProc()
{
  cout << " LightCurveProc::~LightCurveProc() Writing DT-Stars, DT-TimeInfo to file \n" 
       << _outname << endl; 
  POutPersist po(_outname); 
  po << PPFNameTag("stardt") << dts;
  po << PPFNameTag("timdt") << dtti;
}  

void LightCurveProc::ProcessLC(int numet, int nmes, STARINFO *sti, TIMEINFO *tminf, MESUREU *mesu)
{
  double xnt[20];
  xnt[0] = sti->NumEt;
  xnt[1] = sti->FluxRef;
  xnt[2] = sti->XPos;
  xnt[3] = sti->YPos;

  xnt[4] = (double)nmes;

  // We have to check that  fgbad has been filled
  if (fgbad.size() < nmes) 
    for(int kb=fgbad.size(); kb<nmes; kb++) fgbad.push_back(0);

  xnt[5] = 0.; 
  for(int jt=0; jt<nmes; jt++) { /* mesure jt=0 is the reference image */
    if (fgbad[jt] != 0)  continue;
    xnt[5] += mesu[jt].Flux; 
  }                   
  xnt[5] /= nmes;  // mean flux of the star number numet

  xnt[6] = 0.; 
  for(int jt=1; jt<nmes; jt++) {
    if (fgbad[jt] != 0)  continue;
    xnt[6] += (mesu[jt].Flux - xnt[5])*(mesu[jt].Flux - xnt[5]); 
  }
  xnt[6] = sqrt(xnt[6]/nmes); // flux dispersion of the star number numet

  dts.AddRow(xnt);
}

void LightCurveProc::TimeInfo(int nmes, TIMEINFO *tminf)
{
  double xnt[20]; 
  fgbad.clear();
  for(int jt=0; jt<nmes; jt++){
    fgbad.push_back(0);
    if (jt == 0) fgbad[jt] = 1;  // flag mes #0 (ref) as bad 
    if((jt==260)||(jt==647)) fgbad[jt] = 2;  // in-proper mesures
    xnt[0] = tminf[jt].TStart/86400.;
    xnt[1] = tminf[jt].Expose;
    xnt[2] = tminf[jt].Fond;
    xnt[3] = tminf[jt].SigFond;
    xnt[4] = (double)tminf[jt].SigX;
    xnt[5] = tminf[jt].SigY;
    xnt[6] = tminf[jt].Absorption;
    xnt[7] = fgbad[jt];
    dtti.AddRow(xnt);
  }
}