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 |
|
---|
12 | FlagOutOfRange::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 |
|
---|
21 | FlagOutOfRange::~FlagOutOfRange()
|
---|
22 | {
|
---|
23 | }
|
---|
24 |
|
---|
25 | void FlagOutOfRange::Print(::ostream & os)
|
---|
26 | {
|
---|
27 | os<<"FlagOutOfRange::Print Range: ["<<VMin<<","<<VMax<<"] flgBad="<<flgBad<<endl;
|
---|
28 | }
|
---|
29 |
|
---|
30 | void FlagOutOfRange::init()
|
---|
31 | {
|
---|
32 | cout << "FlagOutOfRange::init" << endl;
|
---|
33 | declareInput("DataIn"); // input index 0
|
---|
34 | declareOutput("DataOut"); // output index 0
|
---|
35 | }
|
---|
36 |
|
---|
37 | void FlagOutOfRange::run()
|
---|
38 | {
|
---|
39 | long snb = getMinIn();
|
---|
40 | long sne = getMaxIn();
|
---|
41 |
|
---|
42 | if(snb>sne) {
|
---|
43 | cout<<"FlagOutOfRange::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
|
---|
44 | throw ParmError("FlagOutOfRange::run() - Bad sample interval");
|
---|
45 | }
|
---|
46 | if(!checkInputTOIIndex(0)) {
|
---|
47 | cout<<"FlagOutOfRange::run() - Input TOI (DataIn) not connected! "<<endl;
|
---|
48 | throw ParmError("FlagOutOfRange::run() Input TOI (DataIn) not connected!");
|
---|
49 | }
|
---|
50 | if(!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 | //---------------------------------------------------------
|
---|
56 | try {
|
---|
57 |
|
---|
58 | cout<<"FlagOutOfRange::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
|
---|
59 | uint_4 nbad=0;
|
---|
60 |
|
---|
61 | for(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 | }
|
---|
70 | cout<<"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 |
|
---|
79 | return;
|
---|
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 |
|
---|
97 | FlagGlitch::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 |
|
---|
113 | FlagGlitch::~FlagGlitch()
|
---|
114 | {
|
---|
115 | }
|
---|
116 |
|
---|
117 | void 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 | ////////////////////////////////////////////////////////////////////////
|
---|
126 | void 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 |
|
---|
136 | void FlagGlitch::run()
|
---|
137 | {
|
---|
138 | long snb = getMinIn();
|
---|
139 | long sne = getMaxIn();
|
---|
140 |
|
---|
141 | if(snb>sne) {
|
---|
142 | cout<<"FlagGlitch::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
|
---|
143 | throw ParmError("FlagGlitch::run() - Bad sample interval");
|
---|
144 | }
|
---|
145 | if(!checkInputTOIIndex(0)) {
|
---|
146 | cout<<"FlagGlitch::run() - Input TOI (DataIn) not connected! "<<endl;
|
---|
147 | throw ParmError("FlagGlitch::run() Input TOI (DataIn) not connected!");
|
---|
148 | }
|
---|
149 | if(!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 | //---------------------------------------------------------
|
---|
155 | try {
|
---|
156 |
|
---|
157 | SLinParBuff Slpb0(LSide,BuffUpdate,0.,0.,true);
|
---|
158 | SLinParBuff Slpb1(LSide,BuffUpdate,0.,0.,true);
|
---|
159 | SLinParBuff* Slpb[2] = {&Slpb0,&Slpb1};
|
---|
160 |
|
---|
161 | char* nament_[8] = {"s","bol","sg","sd","slg","sld","ypg","ypd"};
|
---|
162 | r_4 xnt_[8];
|
---|
163 | NTuple nt_(8,nament_);
|
---|
164 | int intuple = mLP/10;
|
---|
165 |
|
---|
166 | if(mLP) cout<<"FlagGlitch::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
|
---|
167 | uint_4 nflag=0;
|
---|
168 | r_8 nsg = NSCut;
|
---|
169 | int_4 lg = LGlitch;
|
---|
170 | int_4 lgs2 = (lg-1)/2;
|
---|
171 | int_4 lm = LSide;
|
---|
172 | int_4 lmmin = lm/2; if(lmmin<2) lmmin=2;
|
---|
173 |
|
---|
174 | // Initialisation
|
---|
175 | if(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 |
|
---|
181 | for(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 | }
|
---|
283 | if(mLP) cout<<"FlagGlitch::run(): Fin de boucle: nflag="<<nflag
|
---|
284 | <<" / "<<sne-snb+1<<endl;
|
---|
285 |
|
---|
286 | if(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 |
|
---|
298 | return;
|
---|
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 |
|
---|
309 | FlagAroundFlag::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 |
|
---|
319 | FlagAroundFlag::~FlagAroundFlag()
|
---|
320 | {
|
---|
321 | }
|
---|
322 |
|
---|
323 | void 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 |
|
---|
331 | void 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 |
|
---|
341 | void FlagAroundFlag::run()
|
---|
342 | {
|
---|
343 | long snb = getMinIn();
|
---|
344 | long sne = getMaxIn();
|
---|
345 |
|
---|
346 | if(snb>sne) {
|
---|
347 | cout<<"FlagAroundFlag::run() - Bad sample interval"<<snb<<" , "<<sne<<endl;
|
---|
348 | throw ParmError("FlagAroundFlag::run() - Bad sample interval");
|
---|
349 | }
|
---|
350 | if(!checkInputTOIIndex(0)) {
|
---|
351 | cout<<"FlagAroundFlag::run() - Input TOI (DataIn) not connected! "<<endl;
|
---|
352 | throw ParmError("FlagAroundFlag::run() Input TOI (DataIn) not connected!");
|
---|
353 | }
|
---|
354 | if(!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 | //---------------------------------------------------------
|
---|
360 | try {
|
---|
361 |
|
---|
362 | cout<<"FlagAroundFlag::run(): SampleNum de "<<snb<<" a "<<sne<<endl;
|
---|
363 | uint_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
|
---|
367 | for(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 | }
|
---|
388 | cout<<"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 |
|
---|
398 | return;
|
---|
399 | }
|
---|