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

Last change on this file since 2311 was 1762, checked in by aubourg, 24 years ago

toujours et encore magique

File size: 11.4 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// Initialisation
175if(LSide>0) for(int_4 sn=snb;sn<snb+LSide && sn<=sne;sn++) {
176 r_8 bolog; uint_8 fgbolog=0;
177 getData(0,sn,bolog,fgbolog);
178 if(!(fgbolog & flgNotLookAt)) Slpb[1]->Push((r_8)sn,bolog);
179}
180
181for(int_4 sn=snb;sn<=sne;sn++) {
182 int_4 lp = (mLPsnb<=sn && sn<=mLPsne) ? mLP%10: 0;
183 r_8 bolog; uint_8 fgbolog=0;
184
185 getData(0,sn,bolog,fgbolog);
186 if(fgbolog & flgNotLookAt) {
187 putData(0,sn,bolog,fgbolog);
188 continue;
189 }
190
191 // i
192 // <---lm---><--lg--><---lm--->
193 // | c=0 | | c=1 |
194 // ideb | | ifin
195 // ifin ideb
196 r_8 bolo; uint_8 fgbolo=0;
197 int_4 ifill[2] = {sn-lgs2-1,sn+lgs2+lm};
198 if(lp>1) printf(">>>> sn=%d v=(%f,%ld) ifill=[%d,%d]\n"
199 ,sn,bolog,fgbolog,ifill[0],ifill[1]);
200
201 if(ifill[0]<snb) {
202 putData(0,sn,bolog,fgbolog);
203 continue;
204 }
205 getData(0,ifill[0],bolo,fgbolo);
206 if(fgbolo&flgNotLookAt) Slpb[0]->Pop();
207 else Slpb[0]->Push((r_8)ifill[0],bolo);
208 if((int_4)Slpb[0]->NPoints()<lmmin) {
209 putData(0,sn,bolog,fgbolog);
210 continue;
211 }
212
213 if(ifill[1]<=sne) {
214 getData(0,ifill[1],bolo,fgbolo);
215 if(fgbolo&flgNotLookAt) Slpb[1]->Pop();
216 else Slpb[1]->Push((r_8)ifill[1],bolo);
217 } else Slpb[1]->Pop();
218 if((int_4)Slpb[1]->NPoints()<lmmin) {
219 putData(0,sn,bolog,fgbolog);
220 continue;
221 }
222
223 // Calcul des valeurs pour decider
224 double s[2],slin[2],ypred[2];
225 for(int c=0;c<2;c++) {
226 r_8 mean,a0,a1;
227 s[c] = Slpb[c]->Compute(mean);
228 slin[c] = Slpb[c]->Compute(a0,a1);
229 ypred[c] = a0+a1*sn;
230 if(lp>1) printf("... c=%d m=%f s=%f ypred=%f slin=%f (%d)\n"
231 ,c,mean,s[c],ypred[c],slin[c],Slpb[c]->NPoints());
232 if(slin[c]<0. || s[c]<0.) {
233 putData(0,sn,bolog,fgbolog);
234 continue;
235 }
236 }
237
238 if(intuple>0) {
239 xnt_[0] = sn; xnt_[1] = bolog;
240 xnt_[2] = s[0]; xnt_[3] = s[1];
241 xnt_[4] = slin[0]; xnt_[5] = slin[1];
242 xnt_[6] = ypred[0]; xnt_[7] = ypred[1];
243 nt_.Fill(xnt_);
244 }
245
246 // Si les predictions ne sont pas compatibles (sig), pas de deglitch
247 r_8 yy = ypred[1]-ypred[0];
248 r_8 sy = nsg*sqrt(s[0]*s[0]+s[1]*s[1]);
249 if(lp>1) printf("... yy=|%f| > sy=%f ???\n",yy,sy);
250 if(fabs(yy) > sy) {
251 putData(0,sn,bolog,fgbolog);
252 continue;
253 }
254
255 // Si la valeur est entre les 2 predictions, pas de deglitch
256 r_8 dp0 = ypred[0];
257 r_8 dp1 = ypred[1];
258 if(ypred[1]<ypred[0]) {dp0 = ypred[1]; dp1 = ypred[0];}
259 if(lp>1) printf("... %f dans [%f,%f] ???\n",bolog,dp0,dp1);
260 if(dp0<bolog && bolog<dp1) {
261 putData(0,sn,bolog,fgbolog);
262 continue;
263 }
264
265 // A + de N sigmas (lin) des 2 predictions
266 dp0 = fabs(bolog-ypred[0]); slin[0] *= nsg;
267 dp1 = fabs(bolog-ypred[1]); slin[1] *= nsg;
268 if(lp>1) printf("... dp0=%f > %f et dp1=%f > %f ???\n"
269 ,dp0,slin[0],dp1,slin[1]);
270 if(dp0 < slin[0] || dp1 < slin[1]) {
271 putData(0,sn,bolog,fgbolog);
272 continue;
273 }
274
275 // C'est un glitch
276 nflag++;
277 fgbolog |= flgGlitch;
278 if(lp>1)
279 printf("......... Glitch en sn=%d ... %d glitches found\n",sn,nflag);
280 putData(0,sn,bolog,fgbolog);
281
282}
283if(mLP) cout<<"FlagGlitch::run(): Fin de boucle: nflag="<<nflag
284 <<" / "<<sne-snb+1<<endl;
285
286if(intuple>0) {
287 char *fileout_ = "flagglitch.ppf";
288 string tag_ = "nfg";
289 POutPersist pos_(fileout_); pos_.PutObject(nt_,tag_);
290}
291
292//---------------------------------------------------------
293} catch (PException & exc) {
294 cout<<"FlagGlitch: Catched Exception "<<(string)typeid(exc).name()
295 <<"\n .... Msg= "<<exc.Msg()<<endl;
296}
297
298return;
299}
300
301
302////////////////////////////////////////////////////////////////////////
303// Flag des sample autour d'un sample deja flaggue :
304// INPUT:
305// lm : longueur a flagguer de part et d'autre du mauvais sample
306// <---lm---><bad_sample><---lm--->
307////////////////////////////////////////////////////////////////////////
308
309FlagAroundFlag::FlagAroundFlag(uint_4 lm,uint_8 flgs)
310{
311 // Mise en forme des arguments
312 LSide = lm;
313
314 SetFlag();
315 SetFlagAroundFlag(flgs);
316 SetLimits();
317}
318
319FlagAroundFlag::~FlagAroundFlag()
320{
321}
322
323void FlagAroundFlag::Print(::ostream & os)
324{
325 os<<"FlagAroundFlag::Print"
326 <<" LSide="<<LSide
327 <<" flgSample="<<flgSample<<" flgAround="<<flgAround<<endl
328 <<" ... Range: ["<<VMin<<","<<VMax<<"]"<<" flgBad="<<flgBad<<endl;
329}
330
331void FlagAroundFlag::init()
332{
333 cout << "FlagAroundFlag::init" << endl;
334 declareInput("DataIn"); // input index 0
335 declareOutput("DataOut"); // output index 0
336
337 int_4 lbuff = 2*LSide +10;
338 if(neededHistory<lbuff) setNeededHistory(lbuff);
339}
340
341void FlagAroundFlag::run()
342{
343long snb = getMinIn();
344long sne = getMaxIn();
345
346if(snb>sne) {
347 cout<<"FlagAroundFlag::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
348 throw ParmError("FlagAroundFlag::run() - Bad sample interval");
349}
350if(!checkInputTOIIndex(0)) {
351 cout<<"FlagAroundFlag::run() - Input TOI (DataIn) not connected! "<<endl;
352 throw ParmError("FlagAroundFlag::run() Input TOI (DataIn) not connected!");
353}
354if(!checkOutputTOIIndex(0)) {
355 cout<<"FlagAroundFlag::run() - Output TOI (DataOut) not connected! "<<endl;
356 throw ParmError("FlagAroundFlag::run() Output TOI (DataOut) not connected!");
357}
358
359//---------------------------------------------------------
360try {
361
362cout<<"FlagAroundFlag::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
363uint_4 nflag=0, nflagged=0, nflagalready=0, nbad=0;
364
365// Un peu lourdingue! On pourrait optimiser!
366// Si on ne tue pas (LSide=0), alors juste un transfert de donnees
367for(int_4 sn=snb;sn<=sne;sn++) {
368 r_8 bolo; uint_8 fgbolo=0;
369 getData(0,sn,bolo,fgbolo);
370 if(VMin<=VMax && (bolo<VMin || bolo>VMax)) { // Pixel out of range
371 nbad++;
372 fgbolo |= flgBad;
373 } else if(fgbolo&flgSample) { // Pixel flaggue flgSample
374 nflag++;
375 } else if(fgbolo&flgAround) { // Pixel deja flagge around?
376 nflagalready++;
377 } else if(LSide>0) { // Pixel a flagguer?
378 int_4 ideb = sn-LSide, ifin = sn+LSide;
379 for(int_4 i=ideb;i<=ifin;i++) { // Y a t-il un "mauvais" pixel autour?
380 if(i<snb || i>sne || i==sn) continue;
381 r_8 b; uint_8 fg=0;
382 getData(0,i,b,fg);
383 if(fg&flgSample) {fgbolo |= flgAround; nflagged++; break;}
384 }
385 }
386 putData(0,sn,bolo,fgbolo);
387}
388cout<<"FlagAroundFlag::run(): Fin de boucle: nflag="<<nflag
389 <<" nflagged="<<nflagged<<" (already="<<nflagalready<<")"<<endl
390 <<" ...... bad="<<nbad<<" / tot="<<sne-snb+1<<endl;
391
392//---------------------------------------------------------
393} catch (PException & exc) {
394 cout<<"FlagAroundFlag: Catched Exception "<<(string)typeid(exc).name()
395 <<"\n .... Msg= "<<exc.Msg()<<endl;
396}
397
398return;
399}
Note: See TracBrowser for help on using the repository browser.