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

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

////////////////////////////////////////////////////////////////////////
// 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("BoloIn");       // input index 0
  declareOutput("BoloOut");     // 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 (BoloIn) not connected! "<<endl;
  throw ParmError("FlagGlitch::run() Input TOI (BoloIn) not connected!");
}
if(!checkOutputTOIIndex(0)) {
  cout<<"FlagGlitch::run() - Output TOI (BoloOut) not connected! "<<endl;
  throw ParmError("FlagGlitch::run() Output TOI (BoloOut) 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;                                                                            
}