source: Sophya/trunk/FrEROS/AnaLC/nbsreadu_base.cc@ 3822

/*
- Test des divers calcul de ligne de base pour la couleur 0
                           (la couleur 1 n'est pas traitee)
*/
#include "sopnamsp.h"
#include "machdefs.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

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

#include "ntuple.h"

void FILL_USER(void);

#define N 5

float *TimeU,*FluxU,*FluxFlU,*eFluxU,*eFluxFlU,*FluxUnoA,*AmpliU;
int_4 npt_;

POutPersist *pos_ = NULL;
NTuple *ntbase = NULL;

/*==========================================================================*/
void UINIT(void)
{
 pos_ = new POutPersist("nbsread_base.ppf");
}

/*===========================================================================*/
#define NXNT 28
void UINITCCD(void)
{
int_4 n;
n = ( nmesure[0] > 0 ) ? nmesure[0] : 1 ;
TimeU     = (float*) malloc((size_t) n*sizeof(float));
FluxU     = (float*) malloc((size_t) n*sizeof(float));
eFluxU    = (float*) malloc((size_t) n*sizeof(float));
FluxFlU   = (float*) malloc((size_t) n*sizeof(float));
eFluxFlU  = (float*) malloc((size_t) n*sizeof(float));
FluxUnoA  = (float*) malloc((size_t) n*sizeof(float));
AmpliU    = (float*) malloc((size_t) n*sizeof(float));

// Create NTuple
const char *nament[NXNT] = {
      "et","xref","flr","dm1","dm21","xp","yp","npt"
     ,"mean","sig","sigi"
     ,"meanB","sigB"
     ,"meanS","sigS","rcS"
     ,"meanSs","sigSs","rcSs"
     ,"meanP","rcP"
     ,"meanPs","rcPs"
     ,"u0","t0","tau","a0","npta"
};
if(ntbase == NULL) delete ntbase;
ntbase = new NTuple(NXNT,nament);

}

/*=========================================================================*/
void UDATCLEAN(int coul)
/* pour tuer une photo selon des criteres utilisateur, mettre date a GRAND2 */
{
int_4 ic = coul-1;
if(nmes[ic]<=0) return;
for(int i=0;i<nmes[ic];i++) {
  if(date[ic][i]<0.)  date[ic][i] = GRAND2;
  if(timeu[ic][i].FgCalib<=0) date[ic][i] = GRAND2;
}
}

/*=========================================================================*/
void UEVT(void)
{
int i,n,rc,nptas;
float mean,sigmaI,sigma,most;
float xnt[NXNT];

/* association rouge-bleu existante et correcte */
if(staru[0].NumEt <= 0 ) return;
if(staru[0].FluxRef<= 0.) return;
for(i=0;i<NXNT;i++) xnt[i] = -9999.;

FILL_USER();

if( iet[0]%50 == 0 )
  printf("UEVT: et=[%7d,%7d] pt nmesure=[%5d] npt=[%5d]\n"
        ,iet[0],iet[1],nmesure[0],npt_);

/* on remplit le ntuple de selection */
xnt[0]=  iet[0];
xnt[1]=  staru[0].XRef;
xnt[2] = staru[0].FluxRef;
xnt[3] = staru[0].DisMin;
xnt[4] = staru[0].DisM2;
xnt[5] = staru[0].XPos;
xnt[6] = staru[0].YPos;
xnt[7] = npt_;

/****************** Filtre Median *******************************/
n = npt_;
rc = FiltMed(FluxU,eFluxU,FluxFlU,eFluxFlU,n);
if(rc<0) for(i=0;i<npt_;i++) {FluxFlU[i]=FluxU[i]; eFluxFlU[i]=eFluxU[i];}

/****************** Moyenne non amplifiee *******************************/
n = npt_;
rc = MeanLine(FluxUnoA,eFluxU,&n,2,4.,&mean,&sigma);
if(rc==0) {
  xnt[8] = mean;
  xnt[9] = sigma;
}

/*********************** Calcul du sigma interne ***************************/
nptas = 0;
n = npt_;
sigmaI = SigmaInt(TimeU,FluxU,eFluxU,&n,4.);
if(sigmaI>0) {
  xnt[10] = sigmaI;
}
if(sigmaI>0 && mean>0) {
  for(i=0;i<npt_;i++) if(eFluxU[i]>0 && (AmpliU[i]-1.)*mean<sigmaI) nptas++;
}

/****************** Calcul de la base sur donnees filtrees ************/
n = npt_; rc=-1;
if(sigmaI>0) rc = BaseLine(FluxFlU,eFluxFlU,&n,4.,sigmaI,&mean,&sigma,&most);
if(rc>=0) {
  xnt[11] = most;
  xnt[12] = sigma;
}

rc = BaseLineS(FluxU,eFluxU,npt_,N,&mean,&sigma);
if(rc>0) {
  xnt[13] = mean;
  xnt[14] = sigma;
  xnt[15] = rc;
}

rc = BaseLineS(FluxU,eFluxU,npt_,-N,&mean,&sigma);
if(rc>0) {
  xnt[16] = mean;
  xnt[17] = sigma;
  xnt[18] = rc;
}

rc = BaseLineP(FluxU,eFluxU,npt_,N,&mean);
if(rc>0) {
  xnt[19] = mean;
  xnt[20] = rc;
}

rc = BaseLineP(FluxU,eFluxU,npt_,-N,&mean);
if(rc>0) {
  xnt[21] = mean;
  xnt[22] = rc;
}

xnt[23] = mc.U0Sim;
xnt[24] = mc.T0Sim;
xnt[25] = mc.TauSim;
xnt[26] = mc.A0Max;
xnt[27] = nptas;
 
ntbase->Fill(xnt);

}

/*=========================================================================*/
void UENDCCD(void)
{
free(TimeU); free(FluxU); free(eFluxU); free(FluxUnoA);
free(FluxFlU); free(eFluxFlU); free(AmpliU);

if(ntbase->NEntry()>0) pos_->PutObject(*ntbase,"ntbase");
if(ntbase!=NULL) delete ntbase; ntbase = NULL;
}

/*=========================================================================*/
void UEND(void)
{
 if(pos_!=NULL) delete pos_; pos_ = NULL;
}

/*=========================================================================*/
void FILL_USER(void)
{
int i,j;
double Flux,FluxB,Xi2,ErrFlux;

npt_=0;
for(i=0;i<nmesure[0];i++) {
  j = indexu[0][i];
  FluxB   = mesu[0][j].FluxB;
  Flux    = mesu[0][j].Flux;
  ErrFlux = mesu[0][j].ErrFlux;
  Xi2     = mesu[0][j].Xi2;
  if( ErrFlux<0. ) ErrFlux *= -1.;
  if( Xi2<=0. || ErrFlux==0. ) continue;
  TimeU[npt_] = date[0][j];
  FluxUnoA[npt_] = Flux;
  if(mc.montecar > 0) Flux *= ampli[0][j];
  AmpliU[npt_] = ampli[0][j];
  FluxU[npt_] = Flux;
  eFluxU[npt_] = ErrFlux;
  npt_++;
}

}