source: Sophya/trunk/FrEROS/AnaLC/nbsreadu_ppf.cc@ 3572

/*
- Creation de NTuple et de courbes de lumieres:
nbsreadu_squel fichier.dataccd -prt 1 \
           -mc 2 -u0 0.01,0.5 -t0 1,-1,-1 -tau 2,10,-1 -tf 2,0,2 -bl 1
*/
#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

#include "fsvcache.h"
#include "fsvst.h"
#include "fmath.h"
#include "nbmath.h"
#include "nbrandom.h"
#include "nbsread.h"
#include "nbgene.h"

#include "ntuple.h"
#include "tvector.h"

POutPersist *pos_ = NULL;
NTuple *ntsti_[2] = {NULL,NULL};

/*================================================================*/
void UINIT(void)
{
 printf("UINIT\n");
 pos_ = new POutPersist("nbsread_ppf.ppf");
}

/*=================== ============================================*/
void UDATCLEAN(int coul)
/*
  Attention ici coul = 1,2
    donc pour tuer une image pour la couleur coul:
    date[coul-1][...] = GRAND2;
*/
{
 int ic = coul - 1;
 cout<<"UDATCLEAN: couleur "<<coul<<", indice tableau "<<ic<<endl;
 if(nmes[ic]<=0) return;

 TVector<r_8> vdate(nmes[ic]); vdate = 0;
 for(int i=0;i<nmes[ic];i++) {
   vdate(i) = date[ic][i];
   if(date[ic][i]<0.)  date[ic][i] = GRAND2;
   if(timeu[ic][i].FgCalib<=0) date[ic][i] = GRAND2;
 }

 char str[16]; sprintf(str,"date%d",ic);
 pos_->PutObject(vdate,str);

}

/*================================================================*/
void UINITCCD(void)
{
 printf("UINITCCD: Debut de l'etude du CCD %d %d\n",numccd,ccdnum);

 // Fill TimeInfoU into NTuple
 const int nxtn = 12;
 const char *nament[nxtn] =
   {"ts","texp","fd","sfd","sx","sy","rho","dx","dy","absor","amass","fgcal"};
 r_8 xnt[nxtn];
 NTuple nttinf(nxtn,nament);
 
 for(int ic=0;ic<NCOULMX;ic++) {
   if(nmesure[ic]<=0) continue;
   NTuple nttinf(nxtn,nament);
   for(int j=0;j<nmesure[ic];j++) {
     int i = indexu[ic][j];
     xnt[0] = timeu[ic][i].TStart/86400.;
     xnt[1] = timeu[ic][i].Expose;
     xnt[2] = timeu[ic][i].Fond;
     xnt[3] = timeu[ic][i].SigFond;
     xnt[4] = timeu[ic][i].SigX;
     xnt[5] = timeu[ic][i].SigY;
     xnt[6] = timeu[ic][i].Rho;
     xnt[7] = timeu[ic][i].DelX;
     xnt[8] = timeu[ic][i].DelY;
     xnt[9] = timeu[ic][i].Absorption;
     xnt[10] = timeu[ic][i].AirMass;
     xnt[11] = timeu[ic][i].FgCalib;
     nttinf.Fill(xnt);
   }
   char str[16]; sprintf(str,"tinf%d",ic);
   pos_->PutObject(nttinf,str);
 }

 // Create NTuple for StarInfo
 const int nsti = 12;
 const char *namentsti[nsti] =
   {"net","xref","fmean","fsig","fref","xp","yp","dm","dm2","dm2r","fgref","nbvois"};
 for(int ic=0;ic<NCOULMX;ic++) {
   cout<<"...Creating NTuple for StarInfo ic="<<ic<<endl;
   if(ntsti_[ic] == NULL) delete ntsti_[ic];
   ntsti_[ic] = new NTuple(nsti,namentsti);
 }

}

/*================================================================*/
void UEVT(void)
{
 printf("****> Et=%d (%d) %d (%d)  fref=%g %g  nmes=%d %d  nmesure=%d %d\n"
       ,iet[0],staru[0].XRef,iet[1],staru[1].XRef,staru[0].FluxRef,staru[1].FluxRef
       ,nmes[0],nmes[1],nmesure[0],nmesure[1]);

 // Fill StarInfo information into NTuple
 for(int ic=0;ic<NCOULMX;ic++) {
   if(iet[ic]<=0) continue;   // pas d'etoile associee
   int n = ntsti_[ic]->NVar();
   if(n<=0) continue;
   r_8 *xnt = new r_8[n];
   xnt[0] = staru[ic].NumEt;
   xnt[1] = staru[ic].XRef;
   xnt[2] = staru[ic].FlxMean;
   xnt[3] = staru[ic].FlxSig;
   xnt[4] = staru[ic].FluxRef;
   xnt[5] = staru[ic].XPos;
   xnt[6] = staru[ic].YPos;
   xnt[7] = staru[ic].DisMin;
   xnt[8] = staru[ic].DisM2;
   xnt[9] = staru[ic].DisM2R;
   xnt[10] = staru[ic].FgRef;
   xnt[11] = staru[ic].NbVois;
   ntsti_[ic]->Fill(xnt);
   delete [] xnt;
 }

 // Fill light curve into NTuple
 if(FgSelSt==iet[0]) {
   cout<<"...Writing light curve for "<<iet[0]<<" "<<iet[1]<<endl;

   const int nxtn = 15;
   const char *nament[nxtn] =
     {"ic",
      "ts","fd","sx","sy","rho",
      "f","fb","xi2","ef","efb","fdl","s9p","pxm",
      "fr"
     };
   r_8 xnt[nxtn];
   NTuple nt(nxtn,nament);

   for(int ic=0;ic<NCOULMX;ic++) {
     if(iet[ic]<=0) continue;   // pas d'etoile associee
     if(nmesure[ic]<=0) continue;
     for(int j=0;j<nmesure[ic];j++) {
       int i = indexu[ic][j];
       xnt[0] = ic;
       xnt[1] = timeu[ic][i].TStart/24.;
       xnt[2] = timeu[ic][i].Fond;
       xnt[3] = timeu[ic][i].SigX;
       xnt[4] = timeu[ic][i].SigY;
       xnt[5] = timeu[ic][i].Rho;
       xnt[6] = mesu[ic][i].Flux;
       xnt[7] = mesu[ic][i].FluxB;
       xnt[8] = mesu[ic][i].Xi2;
       xnt[9] = mesu[ic][i].ErrFlux;
       xnt[10] = mesu[ic][i].ErrFluxB;
       xnt[11] = mesu[ic][i].Fond;
       xnt[12] = mesu[ic][i].S9Pix;
       xnt[13] = mesu[ic][i].PixMax;
       xnt[14] = staru[ic].FluxRef;
       nt.Fill(xnt);
     }
   }
   if(nt.NEntry()>0) {
     char str[64]; sprintf(str,"e%d",iet[0]);
     pos_->PutObject(nt,str);
   }

 }

}

/*================================================================*/
void UENDCCD(void)
{
 printf("UENDCCD: Fin de l'etude du CCD %d %d\n",numccd,ccdnum);
 if(ntsti_[0]->NEntry()>0) pos_->PutObject(*ntsti_[0],"sti0");
 if(ntsti_[1]->NEntry()>0) pos_->PutObject(*ntsti_[1],"sti1");
 if(ntsti_[0]!=NULL) delete ntsti_[0]; ntsti_[0] = NULL;
 if(ntsti_[1]!=NULL) delete ntsti_[1]; ntsti_[1] = NULL;
}

/*================================================================*/
void UEND(void)
{
 printf("UEND: Fin du job\n");
 if(pos_!=NULL) delete pos_; pos_ = NULL;
}


/*

openppf nbsread_ppf.ppf

zone 1 2
n/plot date0.val%n val>0
n/plot date1.val%n val>0

#### timeinfo
zone 1 2

n/plot tinf0.fd%ts ! ! "circlemarker3 red"
n/plot tinf1.fd%ts ! ! "circlemarker3 blue"

n/plot tinf0.sx%ts ! ! "circlemarker3 red"
n/plot tinf1.sx%ts ! ! "circlemarker3 blue"

n/plot tinf0.sy%ts ! ! "circlemarker3 red"
n/plot tinf1.sy%ts ! ! "circlemarker3 blue"

n/plot tinf0.rho%ts ! ! "circlemarker3 red"
n/plot tinf1.rho%ts ! ! "circlemarker3 blue"

n/plot tinf0.dx%ts ! ! "circlemarker3 red"
n/plot tinf1.dx%ts ! ! "circlemarker3 blue"

n/plot tinf0.dy%ts ! ! "circlemarker3 red"
n/plot tinf1.dy%ts ! ! "circlemarker3 blue"

n/plot tinf0.absor%ts ! ! "circlemarker3 red"
n/plot tinf1.absor%ts ! ! "circlemarker3 blue"

n/plot tinf0.amass%ts ! ! "circlemarker3 red"
n/plot tinf1.amass%ts ! ! "circlemarker3 blue"

n/plot tinf0.fgcal%ts ! ! "circlemarker3 red"
n/plot tinf1.fgcal%ts ! ! "circlemarker3 blue"

#### starinfo
zone

n/plot sti0.yp%xp ! ! "nsta red crossmarker5"
n/plot sti1.yp%xp  ! ! "nsta blue circlemarker5 same"

n/plot sti0.log10(fref) fref>0 ! ! "red"
n/plot sti1.log10(fref) fref>0 ! ! "blue keepbin same"

n/plot sti0.dm xref>0 ! ! "red"
n/plot sti1.dm xref>0 ! ! "blue keepbin same"

n/plot sti0.dm2 xref>0 ! ! "red"
n/plot sti1.dm2 xref>0 ! ! "blue keepbin same"

n/plot sti0.nbvois ! ! ! "red"
n/plot sti1.nbvois ! ! ! "blue keepbin same"

n/plot sti0.dm2%dm xref>0 ! "nsta red crossmarker5"
n/plot sti1.dm2%dm  xref>0 ! "nsta blue circlemarker5"

#### light curve
set et 2005

zone 1 2
n/plot e$et.f%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.f%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.ef%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.ef%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.fb%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.fb%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.efb%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.efb%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.(f-fb)/f%ts ic==0&&f>0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.(f-fb)/f%ts ic==1&&f>0 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.fd%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.fd%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.fdl%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.fdl%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.s9p%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.s9p%ts ic==1 ! "nsta connectpoints blue circlemarker3"

n/plot e$et.pxm%ts ic==0 ! "nsta connectpoints red circlemarker3"
n/plot e$et.pxm%ts ic==1 ! "nsta connectpoints blue circlemarker3"

# Distributions
zone

n/plot e$et.fb ic==0
n/plot e$et.f ic==0 ! ! "keepbin same red"
n/plot e$et.fb ic==1
n/plot e$et.f ic==1 ! ! "keepbin same blue"

n/plot e$et.fdl ic==0 ! ! "red"
n/plot e$et.fdl ic==1 ! ! "keepbin same blue"

n/plot e$et.s9p ic==0 ! ! "red"
n/plot e$et.s9p ic==1 ! ! "keepbin same blue"

n/plot e$et.pxm ic==0 ! ! "red"
n/plot e$et.pxm ic==1 ! ! "keepbin same blue"

# Correlations
zone

n/plot e$et.f%fb ic==0 ! "red circlemarker3"
func x 0 1000 10000 "same"
n/plot e$et.f%fb ic==1 ! "blue circlemarker3"
func x 0 1000 10000 "same"

n/plot e$et.s9p%fb ic==0 ! "red circlemarker3"
func x 0 1000 10000 "same"
n/plot e$et.s9p%fb ic==1 ! "blue circlemarker3"
func x 0 1000 10000 "same"

n/plot e$et.pxm%fb ic==0 ! "red circlemarker3"
func x 0 1000 10000 "same"
n/plot e$et.pxm%fb ic==1 ! "blue circlemarker3"
func x 0 1000 10000 "same"
func x 0 1000 10000 "same"

 */