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

Last change on this file since 1737 was 1737, checked in by cmv, 24 years ago

details sur FlagGlitch cmv 07/11/2001

File size: 11.2 KB
Line 
1#include "machdefs.h"
2#include "toimanager.h"
3#include "pexceptions.h"
4#include "ntuple.h"
5#include "ktoibad.h"
6
7
8////////////////////////////////////////////////////////////////////////
9// Flag des sample hors dynamique
10////////////////////////////////////////////////////////////////////////
11
12FlagOutOfRange::FlagOutOfRange(r_8 vmin,r_8 vmax,uint_8 flag)
13 : VMin(vmin), VMax(vmax), flgBad(flag)
14{
15 if(VMax<VMin) {
16 cout<<"FlagOutOfRange::run() - Bad range interval "<<VMin<<" , "<<VMax<<endl;
17 throw ParmError("FlagOutOfRange::FlagOutOfRange() - Bad range interval");
18 }
19}
20
21FlagOutOfRange::~FlagOutOfRange()
22{
23}
24
25void FlagOutOfRange::Print(ostream & os)
26{
27 os<<"FlagOutOfRange::Print Range: ["<<VMin<<","<<VMax<<"] flgBad="<<flgBad<<endl;
28}
29
30void FlagOutOfRange::init()
31{
32 cout << "FlagOutOfRange::init" << endl;
33 declareInput("DataIn"); // input index 0
34 declareOutput("DataOut"); // output index 0
35}
36
37void FlagOutOfRange::run()
38{
39long snb = getMinIn();
40long sne = getMaxIn();
41
42if(snb>sne) {
43 cout<<"FlagOutOfRange::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
44 throw ParmError("FlagOutOfRange::run() - Bad sample interval");
45}
46if(!checkInputTOIIndex(0)) {
47 cout<<"FlagOutOfRange::run() - Input TOI (DataIn) not connected! "<<endl;
48 throw ParmError("FlagOutOfRange::run() Input TOI (DataIn) not connected!");
49}
50if(!checkOutputTOIIndex(0)) {
51 cout<<"FlagOutOfRange::run() - Output TOI (DataOut) not connected! "<<endl;
52 throw ParmError("FlagOutOfRange::run() Output TOI (DataOut) not connected!");
53}
54
55//---------------------------------------------------------
56try {
57
58cout<<"FlagOutOfRange::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
59uint_4 nbad=0;
60
61for(int_4 sn=snb;sn<=sne;sn++) {
62 r_8 bolo; uint_8 fgbolo=0;
63 getData(0,sn,bolo,fgbolo);
64 if(bolo<VMin || bolo>VMax) { // Pixel out of range
65 nbad++;
66 fgbolo |= flgBad;
67 }
68 putData(0,sn,bolo,fgbolo);
69}
70cout<<"FlagOutOfRange::run(): Fin de boucle: nbad="<<nbad
71 <<" / tot="<<sne-snb+1<<endl;
72
73//---------------------------------------------------------
74} catch (PException & exc) {
75 cout<<"FlagOutOfRange: Catched Exception "<<(string)typeid(exc).name()
76 <<"\n .... Msg= "<<exc.Msg()<<endl;
77}
78
79return;
80}
81
82
83////////////////////////////////////////////////////////////////////////
84// Recherche des glitches :
85// INPUT:
86// lg : longeur maxi d'un glitch (impair, si pair lg++)
87// doit etre suffisamment long pour que les glitches soient <lg
88// suffisamment court pour ne pas tuer la galaxie
89// lm : longueur de la fenetre glissante (de chaque cote)
90// <---lm---><--lg--><---lm--->
91// doit etre suffisamment long pour avoir une bonne stabilite
92// des calculs (a0,a1) et suffisamment court pour ne pas
93// tuer la galaxie.
94// nsg : nombre de sigmas pour valider un glitch
95////////////////////////////////////////////////////////////////////////
96
97FlagGlitch::FlagGlitch(uint_4 lg,uint_4 lm,r_8 nsg)
98 : mLP(0), mLPsnb(-1), mLPsne(-1)
99{
100 // Mise en forme des arguments
101 if(nsg<=0) nsg = 3.;
102 if(lg<=0) lg=3; if(lg%2==0) lg++;
103 if(lm<3) lm = 5*lg;
104
105 LGlitch = lg;
106 LSide = lm;
107 NSCut = nsg;
108 DoNotLookAt();
109 SetFlag();
110 SetBuffUpd();
111}
112
113FlagGlitch::~FlagGlitch()
114{
115}
116
117void FlagGlitch::Print(ostream & os)
118{
119 os<<"FlagGlitch::Print"
120 <<" LGlitch="<<LGlitch<<" LSide="<<LSide<<" NSCut="<<NSCut
121 <<" flgNotLookAt="<<flgNotLookAt<<" flgGlitch="<<flgGlitch
122 <<" BuffUpdate="<<BuffUpdate<<endl;
123}
124
125////////////////////////////////////////////////////////////////////////
126void FlagGlitch::init()
127{
128 if(mLP) cout << "FlagGlitch::init" << endl;
129 declareInput("DataIn"); // input index 0
130 declareOutput("DataOut"); // output index 0
131
132 int_4 lbuff = 2*LSide + LGlitch +10;
133 if(neededHistory<lbuff) setNeededHistory(lbuff);
134}
135
136void FlagGlitch::run()
137{
138long snb = getMinIn();
139long sne = getMaxIn();
140
141if(snb>sne) {
142 cout<<"FlagGlitch::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
143 throw ParmError("FlagGlitch::run() - Bad sample interval");
144}
145if(!checkInputTOIIndex(0)) {
146 cout<<"FlagGlitch::run() - Input TOI (DataIn) not connected! "<<endl;
147 throw ParmError("FlagGlitch::run() Input TOI (DataIn) not connected!");
148}
149if(!checkOutputTOIIndex(0)) {
150 cout<<"FlagGlitch::run() - Output TOI (DataOut) not connected! "<<endl;
151 throw ParmError("FlagGlitch::run() Output TOI (DataOut) not connected!");
152}
153
154//---------------------------------------------------------
155try {
156
157SLinParBuff Slpb0(LSide,BuffUpdate,0.,0.,true);
158SLinParBuff Slpb1(LSide,BuffUpdate,0.,0.,true);
159SLinParBuff* Slpb[2] = {&Slpb0,&Slpb1};
160
161char* nament_[8] = {"s","bol","sg","sd","slg","sld","ypg","ypd"};
162r_4 xnt_[8];
163NTuple nt_(8,nament_);
164int intuple = mLP/10;
165
166if(mLP) cout<<"FlagGlitch::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
167uint_4 nflag=0;
168r_8 nsg = NSCut;
169int_4 lg = LGlitch;
170int_4 lgs2 = (lg-1)/2;
171int_4 lm = LSide;
172int_4 lmmin = lm/2; if(lmmin<2) lmmin=2;
173
174
175
176for(int_4 sn=snb;sn<=sne;sn++) {
177 int_4 lp = (mLPsnb<=sn && sn<=mLPsne) ? mLP%10: 0;
178 r_8 bolog; uint_8 fgbolog=0;
179
180 getData(0,sn,bolog,fgbolog);
181 if(fgbolog & flgNotLookAt) {
182 putData(0,sn,bolog,fgbolog);
183 continue;
184 }
185
186 // i
187 // <---lm---><--lg--><---lm--->
188 // | c=0 | | c=1 |
189 // ideb | | ifin
190 // ifin ideb
191 r_8 bolo; uint_8 fgbolo=0;
192 int_4 ifill[2] = {sn-lgs2-1,sn+lgs2+lm};
193 if(lp>1) printf(">>>> sn=%d v=(%f,%ld) ifill=[%d,%d]\n"
194 ,sn,bolog,fgbolog,ifill[0],ifill[1]);
195
196 if(ifill[0]<snb) {
197 putData(0,sn,bolog,fgbolog);
198 continue;
199 }
200 getData(0,ifill[0],bolo,fgbolo);
201 if(fgbolo&flgNotLookAt) Slpb[0]->Pop();
202 else Slpb[0]->Push((r_8)ifill[0],bolo);
203 if((int_4)Slpb[0]->NPoints()<lmmin) {
204 putData(0,sn,bolog,fgbolog);
205 continue;
206 }
207
208 if(ifill[1]<=sne) {
209 getData(0,ifill[1],bolo,fgbolo);
210 if(fgbolo&flgNotLookAt) Slpb[1]->Pop();
211 else Slpb[1]->Push((r_8)ifill[1],bolo);
212 } else Slpb[1]->Pop();
213 if((int_4)Slpb[1]->NPoints()<lmmin) {
214 putData(0,sn,bolog,fgbolog);
215 continue;
216 }
217
218 // Calcul des valeurs pour decider
219 double s[2],slin[2],ypred[2];
220 for(int c=0;c<2;c++) {
221 r_8 mean,a0,a1;
222 s[c] = Slpb[c]->Compute(mean);
223 slin[c] = Slpb[c]->Compute(a0,a1);
224 ypred[c] = a0+a1*sn;
225 if(lp>1) printf("... c=%d m=%f s=%f ypred=%f slin=%f (%d)\n"
226 ,c,mean,s[c],ypred[c],slin[c],Slpb[c]->NPoints());
227 if(slin[c]<0. || s[c]<0.) {
228 putData(0,sn,bolog,fgbolog);
229 continue;
230 }
231 }
232
233 if(intuple>0) {
234 xnt_[0] = sn; xnt_[1] = bolog;
235 xnt_[2] = s[0]; xnt_[3] = s[1];
236 xnt_[4] = slin[0]; xnt_[5] = slin[1];
237 xnt_[6] = ypred[0]; xnt_[7] = ypred[1];
238 nt_.Fill(xnt_);
239 }
240
241 // Si les predictions ne sont pas compatibles (sig), pas de deglitch
242 r_8 yy = ypred[1]-ypred[0];
243 r_8 sy = nsg*sqrt(s[0]*s[0]+s[1]*s[1]);
244 if(lp>1) printf("... yy=|%f| > sy=%f ???\n",yy,sy);
245 if(fabs(yy) > sy) {
246 putData(0,sn,bolog,fgbolog);
247 continue;
248 }
249
250 // Si la valeur est entre les 2 predictions, pas de deglitch
251 r_8 dp0 = ypred[0];
252 r_8 dp1 = ypred[1];
253 if(ypred[1]<ypred[0]) {dp0 = ypred[1]; dp1 = ypred[0];}
254 if(lp>1) printf("... %f dans [%f,%f] ???\n",bolog,dp0,dp1);
255 if(dp0<bolog && bolog<dp1) {
256 putData(0,sn,bolog,fgbolog);
257 continue;
258 }
259
260 // A + de N sigmas (lin) des 2 predictions
261 dp0 = fabs(bolog-ypred[0]); slin[0] *= nsg;
262 dp1 = fabs(bolog-ypred[1]); slin[1] *= nsg;
263 if(lp>1) printf("... dp0=%f > %f et dp1=%f > %f ???\n"
264 ,dp0,slin[0],dp1,slin[1]);
265 if(dp0 < slin[0] || dp1 < slin[1]) {
266 putData(0,sn,bolog,fgbolog);
267 continue;
268 }
269
270 // C'est un glitch
271 nflag++;
272 fgbolog |= flgGlitch;
273 if(lp>1)
274 printf("......... Glitch en sn=%d ... %d glitches found\n",sn,nflag);
275 putData(0,sn,bolog,fgbolog);
276
277}
278if(mLP) cout<<"FlagGlitch::run(): Fin de boucle: nflag="<<nflag
279 <<" / "<<sne-snb+1<<endl;
280
281if(intuple>0) {
282 char *fileout_ = "flagglitch.ppf";
283 string tag_ = "nfg";
284 POutPersist pos_(fileout_); pos_.PutObject(nt_,tag_);
285}
286
287//---------------------------------------------------------
288} catch (PException & exc) {
289 cout<<"FlagGlitch: Catched Exception "<<(string)typeid(exc).name()
290 <<"\n .... Msg= "<<exc.Msg()<<endl;
291}
292
293return;
294}
295
296
297////////////////////////////////////////////////////////////////////////
298// Flag des sample autour d'un sample deja flaggue :
299// INPUT:
300// lm : longueur a flagguer de part et d'autre du mauvais sample
301// <---lm---><bad_sample><---lm--->
302////////////////////////////////////////////////////////////////////////
303
304FlagAroundFlag::FlagAroundFlag(uint_4 lm,uint_8 flgs)
305{
306 // Mise en forme des arguments
307 LSide = lm;
308
309 SetFlag();
310 SetFlagAroundFlag(flgs);
311 SetLimits();
312}
313
314FlagAroundFlag::~FlagAroundFlag()
315{
316}
317
318void FlagAroundFlag::Print(ostream & os)
319{
320 os<<"FlagAroundFlag::Print"
321 <<" LSide="<<LSide
322 <<" flgSample="<<flgSample<<" flgAround="<<flgAround<<endl
323 <<" ... Range: ["<<VMin<<","<<VMax<<"]"<<" flgBad="<<flgBad<<endl;
324}
325
326void FlagAroundFlag::init()
327{
328 cout << "FlagAroundFlag::init" << endl;
329 declareInput("DataIn"); // input index 0
330 declareOutput("DataOut"); // output index 0
331
332 int_4 lbuff = 2*LSide +10;
333 if(neededHistory<lbuff) setNeededHistory(lbuff);
334}
335
336void FlagAroundFlag::run()
337{
338long snb = getMinIn();
339long sne = getMaxIn();
340
341if(snb>sne) {
342 cout<<"FlagAroundFlag::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
343 throw ParmError("FlagAroundFlag::run() - Bad sample interval");
344}
345if(!checkInputTOIIndex(0)) {
346 cout<<"FlagAroundFlag::run() - Input TOI (DataIn) not connected! "<<endl;
347 throw ParmError("FlagAroundFlag::run() Input TOI (DataIn) not connected!");
348}
349if(!checkOutputTOIIndex(0)) {
350 cout<<"FlagAroundFlag::run() - Output TOI (DataOut) not connected! "<<endl;
351 throw ParmError("FlagAroundFlag::run() Output TOI (DataOut) not connected!");
352}
353
354//---------------------------------------------------------
355try {
356
357cout<<"FlagAroundFlag::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
358uint_4 nflag=0, nflagged=0, nflagalready=0, nbad=0;
359
360// Un peu lourdingue! On pourrait optimiser!
361// Si on ne tue pas (LSide=0), alors juste un transfert de donnees
362for(int_4 sn=snb;sn<=sne;sn++) {
363 r_8 bolo; uint_8 fgbolo=0;
364 getData(0,sn,bolo,fgbolo);
365 if(VMin<=VMax && (bolo<VMin || bolo>VMax)) { // Pixel out of range
366 nbad++;
367 fgbolo |= flgBad;
368 } else if(fgbolo&flgSample) { // Pixel flaggue flgSample
369 nflag++;
370 } else if(fgbolo&flgAround) { // Pixel deja flagge around?
371 nflagalready++;
372 } else if(LSide>0) { // Pixel a flagguer?
373 int_4 ideb = sn-LSide, ifin = sn+LSide;
374 for(int_4 i=ideb;i<=ifin;i++) { // Y a t-il un "mauvais" pixel autour?
375 if(i<snb || i>sne || i==sn) continue;
376 r_8 b; uint_8 fg=0;
377 getData(0,i,b,fg);
378 if(fg&flgSample) {fgbolo |= flgAround; nflagged++; break;}
379 }
380 }
381 putData(0,sn,bolo,fgbolo);
382}
383cout<<"FlagAroundFlag::run(): Fin de boucle: nflag="<<nflag
384 <<" nflagged="<<nflagged<<" (already="<<nflagalready<<")"<<endl
385 <<" ...... bad="<<nbad<<" / tot="<<sne-snb+1<<endl;
386
387//---------------------------------------------------------
388} catch (PException & exc) {
389 cout<<"FlagAroundFlag: Catched Exception "<<(string)typeid(exc).name()
390 <<"\n .... Msg= "<<exc.Msg()<<endl;
391}
392
393return;
394}
Note: See TracBrowser for help on using the repository browser.