source: Sophya/trunk/ArchTOIPipe/ProcWSophya/ktoibad.cc@ 1736

#include "machdefs.h"
#include "toimanager.h"
#include "pexceptions.h"
#include "ntuple.h"
#include "ktoibad.h"


////////////////////////////////////////////////////////////////////////
// Flag des sample hors dynamique
////////////////////////////////////////////////////////////////////////

FlagOutOfRange::FlagOutOfRange(r_8 vmin,r_8 vmax,uint_8 flag)
  : VMin(vmin), VMax(vmax), flgBad(flag)
{
 if(VMax<VMin) {
   cout<<"FlagOutOfRange::run() - Bad range interval "<<VMin<<" , "<<VMax<<endl;
   throw ParmError("FlagOutOfRange::FlagOutOfRange() - Bad range interval");
 }
}

FlagOutOfRange::~FlagOutOfRange()
{
}

void FlagOutOfRange::Print(ostream & os)
{
 os<<"FlagOutOfRange::Print Range: ["<<VMin<<","<<VMax<<"] flgBad="<<flgBad<<endl;
}

void FlagOutOfRange::init()
{
  cout << "FlagOutOfRange::init" << endl;
  declareInput("DataIn");       // input index 0
  declareOutput("DataOut");     // output index 0
}

void FlagOutOfRange::run()
{
long snb = getMinIn();
long sne = getMaxIn();

if(snb>sne) {
  cout<<"FlagOutOfRange::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
  throw ParmError("FlagOutOfRange::run() - Bad sample interval");
}
if(!checkInputTOIIndex(0)) {
  cout<<"FlagOutOfRange::run() - Input TOI (DataIn) not connected! "<<endl;
  throw ParmError("FlagOutOfRange::run() Input TOI (DataIn) not connected!");
}
if(!checkOutputTOIIndex(0)) {
  cout<<"FlagOutOfRange::run() - Output TOI (DataOut) not connected! "<<endl;
  throw ParmError("FlagOutOfRange::run() Output TOI (DataOut) not connected!");
}

//---------------------------------------------------------
try {

cout<<"FlagOutOfRange::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
uint_4 nbad=0;

for(int_4 sn=snb;sn<=sne;sn++) {
  r_8 bolo; uint_8 fgbolo=0;
  getData(0,sn,bolo,fgbolo);
  if(bolo<VMin || bolo>VMax) { // Pixel out of range
    nbad++;
    fgbolo |= flgBad;
  }
  putData(0,sn,bolo,fgbolo);
}
cout<<"FlagOutOfRange::run(): Fin de boucle: nbad="<<nbad
    <<" / tot="<<sne-snb+1<<endl;

//---------------------------------------------------------
} catch (PException & exc) {
  cout<<"FlagOutOfRange: Catched Exception "<<(string)typeid(exc).name()
      <<"\n .... Msg= "<<exc.Msg()<<endl;
}

return;                                                                            
}


////////////////////////////////////////////////////////////////////////
// Recherche des glitches :
// INPUT:
// lg : longeur maxi d'un glitch (impair, si pair lg++)
//      doit etre suffisamment long pour que les glitches soient <lg
//      suffisamment court pour ne pas tuer la galaxie
// lm : longueur de la fenetre glissante (de chaque cote)
//      <---lm---><--lg--><---lm--->
//      doit etre suffisamment long pour avoir une bonne stabilite
//      des calculs (a0,a1) et suffisamment court pour ne pas
//      tuer la galaxie.
// nsg : nombre de sigmas pour valider un glitch
////////////////////////////////////////////////////////////////////////

FlagGlitch::FlagGlitch(uint_4 lg,uint_4 lm,r_8 nsg)
  : mLP(0), mLPsnb(-1), mLPsne(-1)
{
 // Mise en forme des arguments
 if(nsg<=0) nsg = 3.;
 if(lg<=0) lg=3; if(lg%2==0) lg++;
 if(lm<3) lm = 5*lg;

 LGlitch = lg;
 LSide =  lm;
 NSCut = nsg;
 Slpb[0] = Slpb[1] = NULL;
 DoNotLookAt();
 SetFlag();
 SetBuffUpd();
}

FlagGlitch::~FlagGlitch()
{
 for(int i=0;i<2;i++) {if(Slpb[i]) delete Slpb[i]; Slpb[i]=NULL;}
}

void FlagGlitch::Print(ostream & os)
{
 os<<"FlagGlitch::Print"
   <<" LGlitch="<<LGlitch<<" LSide="<<LSide<<" NSCut="<<NSCut
   <<" flgNotLookAt="<<flgNotLookAt<<" flgGlitch="<<flgGlitch
   <<" BuffUpdate="<<BuffUpdate<<endl;
}

////////////////////////////////////////////////////////////////////////
void FlagGlitch::init()
{
  if(mLP) cout << "FlagGlitch::init" << endl;
  declareInput("DataIn");       // input index 0
  declareOutput("DataOut");     // output index 0

  for(int i=0;i<2;i++)
    Slpb[i] = new SLinParBuff(LSide,BuffUpdate,0.,0.,true);

  int_4 lbuff = 2*LSide + LGlitch +10;
  if(neededHistory<lbuff) setNeededHistory(lbuff);
}

void FlagGlitch::run()
{
long snb = getMinIn();
long sne = getMaxIn();

if(snb>sne) {
  cout<<"FlagGlitch::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
  throw ParmError("FlagGlitch::run() - Bad sample interval");
}
if(!checkInputTOIIndex(0)) {
  cout<<"FlagGlitch::run() - Input TOI (DataIn) not connected! "<<endl;
  throw ParmError("FlagGlitch::run() Input TOI (DataIn) not connected!");
}
if(!checkOutputTOIIndex(0)) {
  cout<<"FlagGlitch::run() - Output TOI (DataOut) not connected! "<<endl;
  throw ParmError("FlagGlitch::run() Output TOI (DataOut) not connected!");
}

//---------------------------------------------------------
try {

char* nament_[8] = {"s","bol","sg","sd","slg","sld","ypg","ypd"};
r_4 xnt_[8];
NTuple nt_(8,nament_);
int intuple = mLP/10;

if(mLP) cout<<"FlagGlitch::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
uint_4 nflag=0;
r_8 nsg = NSCut;
int_4 lg = LGlitch;
int_4 lgs2 = (lg-1)/2;
int_4 lm = LSide;
int_4 lmmin = lm/2; if(lmmin<2) lmmin=2;

for(int_4 sn=snb;sn<=sne;sn++) {
  int_4 lp = (mLPsnb<=sn && sn<=mLPsne) ? mLP%10: 0;
  r_8 bolog; uint_8 fgbolog=0;

  getData(0,sn,bolog,fgbolog);
  if(fgbolog & flgNotLookAt) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  //                 i
  //    <---lm---><--lg--><---lm--->
  //    |   c=0  |        |   c=1  |
  // ideb        |        |        ifin
  //             ifin     ideb
  r_8 bolo; uint_8 fgbolo=0;
  int_4 ifill[2] = {sn-lgs2-1,sn+lgs2+lm};
  if(lp>1) printf(">>>> sn=%d v=(%f,%ld)  ifill=[%d,%d]\n"
                 ,sn,bolog,fgbolog,ifill[0],ifill[1]);

  if(ifill[0]<snb) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }
  getData(0,ifill[0],bolo,fgbolo);
  if(fgbolo&flgNotLookAt) Slpb[0]->Pop();
    else Slpb[0]->Push((r_8)ifill[0],bolo);
  if((int_4)Slpb[0]->NPoints()<lmmin) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  if(ifill[1]<=sne) {
    getData(0,ifill[1],bolo,fgbolo);
    if(fgbolo&flgNotLookAt) Slpb[1]->Pop();
      else Slpb[1]->Push((r_8)ifill[1],bolo);
  } else Slpb[1]->Pop();
  if((int_4)Slpb[1]->NPoints()<lmmin) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  // Calcul des valeurs pour decider
  double s[2],slin[2],ypred[2];
  for(int c=0;c<2;c++) {
    r_8 mean,a0,a1;
    s[c] = Slpb[c]->Compute(mean);
    slin[c] = Slpb[c]->Compute(a0,a1);
    ypred[c] = a0+a1*sn;
    if(lp>1) printf("... c=%d m=%f  s=%f   ypred=%f slin=%f (%d)\n"
                   ,c,mean,s[c],ypred[c],slin[c],Slpb[c]->NPoints());
    if(slin[c]<0. || s[c]<0.) {
      putData(0,sn,bolog,fgbolog);
      continue;
    }
  }

  if(intuple>0) {
    xnt_[0] = sn;  xnt_[1] = bolog;
    xnt_[2] = s[0];     xnt_[3] = s[1];
    xnt_[4] = slin[0];  xnt_[5] = slin[1];
    xnt_[6] = ypred[0]; xnt_[7] = ypred[1];
    nt_.Fill(xnt_);
  }

  // Si les predictions ne sont pas compatibles (sig), pas de deglitch
  r_8 yy = ypred[1]-ypred[0];
  r_8 sy = nsg*sqrt(s[0]*s[0]+s[1]*s[1]);
  if(lp>1) printf("... yy=|%f| > sy=%f ???\n",yy,sy);
  if(fabs(yy) > sy) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  // Si la valeur est entre les 2 predictions, pas de deglitch
  r_8 dp0 = ypred[0];
  r_8 dp1 = ypred[1];
  if(ypred[1]<ypred[0]) {dp0 = ypred[1]; dp1 = ypred[0];}
  if(lp>1) printf("... %f dans [%f,%f] ???\n",bolog,dp0,dp1);
  if(dp0<bolog && bolog<dp1) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  // A + de N sigmas (lin) des 2 predictions
  dp0 = fabs(bolog-ypred[0]); slin[0] *= nsg;
  dp1 = fabs(bolog-ypred[1]); slin[1] *= nsg;
  if(lp>1) printf("... dp0=%f > %f  et   dp1=%f > %f ???\n"
                 ,dp0,slin[0],dp1,slin[1]);
  if(dp0 < slin[0] || dp1 < slin[1]) {
    putData(0,sn,bolog,fgbolog);
    continue;
  }

  // C'est un glitch
  nflag++;
  fgbolog |= flgGlitch;
  if(lp>1)
    printf("......... Glitch en sn=%d ...  %d glitches found\n",sn,nflag);
  putData(0,sn,bolog,fgbolog);

}
if(mLP) cout<<"FlagGlitch::run(): Fin de boucle: nflag="<<nflag
            <<" / "<<sne-snb+1<<endl;

if(intuple>0)  {
 char *fileout_ = "flagglitch.ppf";
 string tag_ = "nfg";
 POutPersist pos_(fileout_); pos_.PutObject(nt_,tag_);
}

//---------------------------------------------------------
} catch (PException & exc) {
  cout<<"FlagGlitch: Catched Exception "<<(string)typeid(exc).name()
      <<"\n .... Msg= "<<exc.Msg()<<endl;
}

return;                                                                            
}


////////////////////////////////////////////////////////////////////////
// Flag des sample autour d'un sample deja flaggue :
// INPUT:
// lm : longueur a flagguer de part et d'autre du mauvais sample
//      <---lm---><bad_sample><---lm--->
////////////////////////////////////////////////////////////////////////

FlagAroundFlag::FlagAroundFlag(uint_4 lm,uint_8 flgs)
{
 // Mise en forme des arguments
 LSide = lm;

 SetFlag();
 SetFlagAroundFlag(flgs);
 SetLimits();
}

FlagAroundFlag::~FlagAroundFlag()
{
}

void FlagAroundFlag::Print(ostream & os)
{
 os<<"FlagAroundFlag::Print"
   <<" LSide="<<LSide
   <<" flgSample="<<flgSample<<" flgAround="<<flgAround<<endl
   <<" ... Range: ["<<VMin<<","<<VMax<<"]"<<" flgBad="<<flgBad<<endl;
}

void FlagAroundFlag::init()
{
  cout << "FlagAroundFlag::init" << endl;
  declareInput("DataIn");       // input index 0
  declareOutput("DataOut");     // output index 0

  int_4 lbuff = 2*LSide +10;
  if(neededHistory<lbuff) setNeededHistory(lbuff);
}

void FlagAroundFlag::run()
{
long snb = getMinIn();
long sne = getMaxIn();

if(snb>sne) {
  cout<<"FlagAroundFlag::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
  throw ParmError("FlagAroundFlag::run() - Bad sample interval");
}
if(!checkInputTOIIndex(0)) {
  cout<<"FlagAroundFlag::run() - Input TOI (DataIn) not connected! "<<endl;
  throw ParmError("FlagAroundFlag::run() Input TOI (DataIn) not connected!");
}
if(!checkOutputTOIIndex(0)) {
  cout<<"FlagAroundFlag::run() - Output TOI (DataOut) not connected! "<<endl;
  throw ParmError("FlagAroundFlag::run() Output TOI (DataOut) not connected!");
}

//---------------------------------------------------------
try {

cout<<"FlagAroundFlag::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
uint_4 nflag=0, nflagged=0, nflagalready=0, nbad=0;

// Un peu lourdingue! On pourrait optimiser!
// Si on ne tue pas (LSide=0), alors juste un transfert de donnees
for(int_4 sn=snb;sn<=sne;sn++) {
  r_8 bolo; uint_8 fgbolo=0;
  getData(0,sn,bolo,fgbolo);
  if(VMin<=VMax && (bolo<VMin || bolo>VMax)) { // Pixel out of range
    nbad++;
    fgbolo |= flgBad;
  } else if(fgbolo&flgSample) { // Pixel flaggue flgSample
    nflag++;
  } else if(fgbolo&flgAround) { // Pixel deja flagge around?
    nflagalready++;
  } else if(LSide>0) { // Pixel a flagguer?
    int_4 ideb = sn-LSide, ifin = sn+LSide;
    for(int_4 i=ideb;i<=ifin;i++) { // Y a t-il un "mauvais" pixel autour?
      if(i<snb || i>sne || i==sn) continue;
      r_8 b; uint_8 fg=0;
      getData(0,i,b,fg);
      if(fg&flgSample) {fgbolo |= flgAround; nflagged++; break;}
    }
  }
  putData(0,sn,bolo,fgbolo);
}
cout<<"FlagAroundFlag::run(): Fin de boucle: nflag="<<nflag
    <<" nflagged="<<nflagged<<" (already="<<nflagalready<<")"<<endl
    <<" ......      bad="<<nbad<<" / tot="<<sne-snb+1<<endl;

//---------------------------------------------------------
} catch (PException & exc) {
  cout<<"FlagAroundFlag: Catched Exception "<<(string)typeid(exc).name()
      <<"\n .... Msg= "<<exc.Msg()<<endl;
}

return;                                                                            
}