source: Sophya/trunk/ArchTOIPipe/ProcWSophya/simoffset.cc@ 2025

Last change on this file since 2025 was 2014, checked in by ansari, 23 years ago

Amelioration de calcul d'offset (possibilite de coupure sur coordonnees
galactique)
Nouveau processeur nettoyeur (SimpleCleaner)
Ajout programme de traitement aksj02.cc , donnees Archeops KS1/KS3

Reza 28/5/2002

File size: 10.9 KB
Line 
1// ArchTOIPipe (C) CEA/DAPNIA/SPP IN2P3/LAL
2// Eric Aubourg
3// Christophe Magneville
4// Reza Ansari
5#include "config.h"
6
7#include "array.h"
8#include "simoffset.h"
9#include <math.h>
10#include "toimanager.h"
11#include "pexceptions.h"
12#include "ctimer.h"
13#include "xntuple.h"
14
15#include "flagtoidef.h"
16
17SimpleOffsetEstimator::SimpleOffsetEstimator(int mwsz, int nptfit, int degpol)
18 : poly((degpol > 1)?degpol:1)
19{
20 SetParams(mwsz, nptfit, degpol);
21 totnscount = 0;
22 totnbblock = 0;
23 bmincut = bmaxcut = -9999.;
24 bgalcut = false;
25 ns_flgcut = ns_bgalcut = 0;
26 SavePolyNTuple();
27}
28
29SimpleOffsetEstimator::~SimpleOffsetEstimator()
30{
31}
32
33void SimpleOffsetEstimator::SetParams(int mwsz, int nptfit, int degpol)
34{
35 mWSz = (mwsz > 8) ? mwsz : 8;
36 nPtFit = (nptfit > degpol+2) ? nptfit : degpol+2;
37 if (nPtFit%2 == 0) nPtFit++;
38 degPol = (degpol > 1)?degpol:1;
39
40}
41
42void SimpleOffsetEstimator::SetBGalCut(double bmin, double bmax)
43{
44 bgalcut = true;
45 bmincut = bmin;
46 bmaxcut = bmax;
47}
48
49void SimpleOffsetEstimator::PrintStatus(::ostream & os)
50{
51 os << "\n ------------------------------------------------------ \n"
52 << " SimpleOffsetEstimator::PrintStatus() - MeanWSize= " << mWSz << " NPtFit="
53 << nPtFit << " DegPoly=" << degPol << " poly.Degre()=" << poly.Degre() << endl;
54 if (bgalcut)
55 os << " bGalCut = Yes , bGalMin = " << bmincut << " bGalMax= " << bmaxcut << endl;
56 else os << " bGalCut = No " << endl;
57 TOIProcessor::PrintStatus(os);
58 os << " ProcessedSampleCount=" << ProcessedSampleCount()
59 << " NS_FlgCut= " << ns_flgcut << " NS_BGalCut= " << ns_bgalcut << endl;
60 os << " ------------------------------------------------------ " << endl;
61}
62
63void SimpleOffsetEstimator::init()
64{
65 cout << "SimpleOffsetEstimator::init" << endl;
66 declareInput("in");
67 declareInput("bgal");
68 declareOutput("offset");
69 declareOutput("out");
70 declareOutput("incopie");
71 declareOutput("bgalcopie");
72 // declareOutput("mean_y");
73 // declareOutput("sig_y");
74 // declareOutput("mean_x");
75 name = "SimpleOffsetEstimator";
76}
77
78void SimpleOffsetEstimator::run()
79{
80 int snb = getMinIn();
81 int sne = getMaxIn();
82
83 bool fgoffset = checkOutputTOIIndex(0);
84 bool fgout = checkOutputTOIIndex(1);
85 bool fgincopie = checkOutputTOIIndex(2);
86 bool fgbgalcopie = checkOutputTOIIndex(3);
87 // bool fgmeany = checkOutputTOIIndex(4);
88 // bool fgsigy = checkOutputTOIIndex(5);
89 // bool fgmeanx = checkOutputTOIIndex(6);
90
91 if (!checkInputTOIIndex(0)) {
92 cerr << " SimpleOffsetEstimator::run() - Input TOI (in) not connected! "
93 << endl;
94 throw ParmError("SimpleOffsetEstimator::run() Input TOI (in) not connected!");
95 }
96 if (bgalcut && !checkInputTOIIndex(1)) {
97 cerr << " SimpleOffsetEstimator::run() - Input TOI bgal not connected! "
98 << endl;
99 throw ParmError("SimpleOffsetEstimator::run() Input TOI bgal not connected!");
100 }
101 if (!fgoffset && !fgout) {
102 cerr << " SimpleOffsetEstimator::run() - No Output TOI (offset/in-offset) connected! "
103 << endl;
104 throw ParmError(" SimpleOffsetEstimator::run() No output TOI (offset/in-offset) connected!");
105 }
106
107 cout << " SimpleOffsetEstimator::run() SNRange=" << snb << " - " << sne << endl;
108
109 // NTuple pour sauvegarde des coeff de poly
110 char * nomsnt[] = {"sncur", "sn0", "meanx", "meany", "sigy",
111 "a0", "a1", "a2", "ycur", "nok", "nbblkok"};
112 XNTuple xntp(0, 11, 0, 0, nomsnt);
113
114 try {
115
116 // Vecteurs pour les donnees et les sorties
117 int wsize = mWSz;
118
119 int hisblk = nPtFit/2+1;
120 Vector vinhist(wsize*hisblk);
121 TVector<uint_8> vfghist(wsize*hisblk);
122 Vector voff(wsize);
123 Vector vout(wsize);
124 Vector vinc(wsize);
125 TVector<uint_8> vfgc(wsize);
126 Vector bgal(wsize);
127
128
129 // Pour le fit
130 Vector errCoef(degPol+1);
131
132 Vector X(nPtFit+1);
133 Vector X0(nPtFit+1);
134 Vector Y(nPtFit+1);
135 Vector YErr(nPtFit+1);
136
137 // Variables diverses
138 int k,kb,j,klast;
139 klast = snb-1;
140 totnbblock = 0;
141
142
143 int nbblkok = 0;
144
145 bool fginiXYdone = false;
146
147 double sn0 = 0.;
148 double nok = 0.;
149 double mean = 0.;
150 double sig = 0.;
151 double meanx = 0.;
152 double sn_last = 0.;
153 double y_last = 0.;
154
155 double last_meanok = 0.;
156 double last_sigok = 1.;
157 double meanx_forlast = 0.;
158 bool fg_last_meansig = false;
159
160 // Boucle sur les sampleNum
161 // 1ere partie, on traite par paquets de wsize
162 int nblk = (sne-snb+1)/wsize;
163 for(kb=0; kb<nblk; kb++) {
164 k = kb*wsize+snb;
165 // for(k=snb;k<=sne-wsize+1;k+=wsize) {
166 // Lecture d'un bloc de donnees
167 Vector vin( vinhist(Range((kb%hisblk)*wsize, 0, wsize) ) );
168 TVector<uint_8> vfg( vfghist(Range((kb%hisblk)*wsize, 0, wsize) ) );
169
170 getData(0, k, wsize, vin.Data(), vfg.Data());
171 if (bgalcut) {
172 getData(1, k, wsize, bgal.Data());
173 if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
174 }
175
176 fg_last_meansig = false;
177
178 if (kb == 0) {
179 last_meanok = vin.Sum()/wsize;
180 last_sigok = vin.SumX2()/wsize - last_meanok*last_meanok;
181 }
182
183 // Calcul moyenne et sigma du bloc
184 nok = 0.; meanx = 0.;
185 mean = 0.; sig = 0.;
186 meanx_forlast = 0.;
187 for(j=0; j<wsize; j++) {
188 meanx_forlast += k+j;
189 if ( vfg(j) ) { ns_flgcut++; continue; }
190 if (bgalcut && (bgal(j) > bmincut) && (bgal(j) < bmaxcut) )
191 { ns_bgalcut++; continue; }
192 mean += vin(j);
193 sig += vin(j)*vin(j);
194 meanx += k+j;
195 nok++;
196 }
197
198
199 if (!fginiXYdone) {
200 if (nok > 0.5) {
201 mean /= nok;
202 meanx /= nok;
203 sig = sig/nok-mean*mean;
204 }
205 X = RegularSequence((kb+0.5)*wsize, (double)wsize);
206 Y = mean;
207 YErr = (nok > 0.5) ? sqrt(mean) : 1.;
208 fginiXYdone = true;
209 }
210
211 if (nok > 3.5) {
212 mean /= nok;
213 meanx /= nok;
214 sig = sig/nok-mean*mean;
215 if (sig < 1.e-10*mean) sig = 1.e-10*mean;
216 if (sig < 1.e-39) sig = 1.e-39;
217 int kk = nbblkok%nPtFit;
218 nbblkok++;
219 Y(kk) = mean;
220 YErr(kk) = sig;
221 X(kk) = meanx;
222 last_meanok = mean;
223 last_sigok = sig;
224 }
225 else {
226 int kk = nbblkok%nPtFit;
227 Y(kk) = last_meanok;
228 YErr(kk) = last_sigok*10.;
229 X(kk) = meanx_forlast/wsize;
230 fg_last_meansig = true;
231 }
232
233 if (((nok>3.5) || fg_last_meansig) && (nbblkok >= nPtFit) ) {
234 // On force le fit a garder une certaine continuite
235 Y(nPtFit) = y_last;
236 double smin, smax;
237 YErr(Range(0,0,nPtFit)).MinMax(smin, smax);
238 if (smax < 1.e-39) smax = 1.e-39;
239 YErr(nPtFit) = smax/10.;
240 X(nPtFit) = sn_last;
241 sn0 = (double)(k-nPtFit*wsize*0.5);
242 X0 = X;
243 X0 -= sn0;
244 Vector polsave(degPol);
245 for(int jj=0; jj<=poly.Degre(); jj++) polsave(jj) = poly[jj];
246 try {
247 poly.Fit(X0,Y,YErr,degPol,errCoef);
248 }
249 catch (CaughtSignalExc& excsig) {
250 cout << " -- simoffset.cc/ catched CaughtSignalExc - Msg= "
251 << excsig.Msg() << " kb=" << kb << " k=" << k << endl;
252 cout << " X0= " << X0 << " Y=" << Y << " YErr=" << YErr << endl;
253 for(int jj=0; jj<=poly.Degre(); jj++) poly[jj] = polsave(jj);
254 }
255 }
256 else {
257 if (nbblkok < 2) {
258 sn0 = k+1;
259 poly[0] = mean;
260 for(int jj=1; jj<=poly.Degre(); jj++) poly[jj] = 0.;
261 }
262 else if (nbblkok < nPtFit) {
263 poly[0] = 0.5*(Y(nbblkok-1)+Y(0));
264 poly[1] = (Y(nbblkok-1) - Y(0))/(X(nbblkok-1)-X(0));
265 sn0 = 0.5*(X(nbblkok-1)+X(0));
266 for(int jj=2; jj<=poly.Degre(); jj++) poly[jj] = 0.;
267 }
268 sn_last = sn0;
269 y_last = poly(sn_last-sn0);
270 }
271 if (ntpoly) { // Remplissage du XNTuple de controle
272 float xnt[12];
273 xnt[0] = k;
274 xnt[1] = sn0;
275 xnt[2] = meanx;
276 xnt[3] = mean;
277 xnt[4] = sig;
278 xnt[5] = poly[0];
279 xnt[6] = poly[1];
280 xnt[7] = poly[2];
281 xnt[8] = poly(k-sn0);
282 xnt[9] = nok;
283 xnt[10] = nbblkok;
284 xntp.Fill(NULL, xnt, NULL, NULL);
285 }
286
287 /*
288 if (nbblkok < 8) {
289 cout << "------ DBG-X " << nbblkok << "," << nok << " degre=" << poly.Degre() << endl;
290 cout << "DBG-A X0=" << X0 << endl;
291 cout << "DBG-A Y=" << Y << endl;
292 cout << "DBG-A YErr=" << YErr << endl;
293 cout << "DBG-A poly= " << poly << endl;
294 }
295 */
296
297 if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
298
299 if (kb < nPtFit/2) continue;
300 Vector vinc;
301 TVector<uint_8> vfgc;
302 int kbs = kb-nPtFit/2;
303 k = kbs*wsize+snb;
304 if (kb == nblk-1) {
305 int wszt = wsize*hisblk;
306 voff.ReSize(wszt);
307 vout.ReSize(wszt);
308 vinc.ReSize(wszt);
309 vfgc.ReSize(wszt);
310 bgal.ReSize(wszt);
311 int jb = 0;
312 for(int kbsc=kbs; kbsc<nblk; kbsc++) {
313 vinc(Range(jb*wsize, 0, wsize)) =
314 vinhist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
315 vfgc(Range(jb*wsize, 0, wsize)) =
316 vfghist(Range((kbsc%hisblk)*wsize, 0, wsize) ) ;
317 jb++;
318 }
319 wsize = wszt;
320 }
321 else {
322 vinc = vinhist(Range((kbs%hisblk)*wsize, 0, wsize) ) ;
323 vfgc = vfghist(Range((kbs%hisblk)*wsize, 0, wsize) ) ;
324 }
325
326 // Calcul des valeurs d'offset en sortie
327 for(j=0; j<wsize; j++)
328 voff(j) = poly(k+j-sn0);
329 sn_last = k+wsize;
330 y_last = poly(sn_last-sn0);
331
332 if (fgoffset) putData(0, k, wsize, voff.Data());
333 if (fgout) {
334 vinc -= voff;
335 putData(1, k, wsize, vinc.Data(), vfgc.Data());
336 }
337 /*
338 if (fgmeany) {
339 vout = mean;
340 putData(4, k, wsize, vout.Data());
341 }
342 if (fgsigy) {
343 vout = sig;
344 putData(5, k, wsize, vout.Data());
345 }
346
347 if (fgmeanx) {
348 vout = meanx;
349 putData(6, k, wsize, vout.Data());
350 }
351 */
352
353 klast+=wsize;
354 totnscount+=wsize;
355 totnbblock++;
356
357 } // Fin boucle sur les samples, par pas de wsize
358
359 // 3eme partie, on traite la fin du bloc d'echantillons si necessaire
360 if (klast < sne) {
361 wsize = sne-klast;
362 Vector vin(wsize);
363 voff.ReSize(wsize);
364 vout.ReSize(wsize);
365 TVector<uint_8> vfg(wsize);
366 k = klast+1;
367 getData(0, k, wsize, vin.Data(), vfg.Data());
368
369 if (bgalcut) {
370 getData(1, k, wsize, bgal.Data());
371 if (fgbgalcopie) putData(3, k, wsize, bgal.Data());
372 }
373
374 for(j=0; j<wsize; j++)
375 voff(j) = poly(k+j-sn0);
376 if (fgoffset) putData(0, k, wsize, voff.Data());
377 if (fgincopie) putData(2, k, wsize, vin.Data(), vfg.Data());
378 if (fgout) {
379 vin -= voff;
380 putData(1, k, wsize, vin.Data(), vfg.Data());
381 }
382
383 /*
384 if (fgmeany) {
385 vout = mean;
386 putData(4, k, wsize, vout.Data());
387 }
388 if (fgsigy) {
389 vout = sig;
390 putData(5, k, wsize, vout.Data());
391 }
392
393 if (fgmeanx) {
394 vout = meanx;
395 putData(6, k, wsize, vout.Data());
396 }
397 */
398
399 klast+=wsize;
400 totnscount+=wsize;
401 totnbblock++;
402 }
403
404 cout << " SimpleOffsetEstimator::run() - End of processing "
405 << " TotNbBlocks= " << totnbblock << " ProcSamples=" << totnscount << endl;
406
407 } // Bloc try
408
409 catch (PException & exc) {
410 cerr << "SimpleOffsetEstimator::run() Catched Exception " << (string)typeid(exc).name()
411 << "\n .... Msg= " << exc.Msg() << endl;
412 }
413
414 if (ntpoly) {
415 if (ntpolyname.length() < 1) ntpolyname = "simoffset.ppf";
416 POutPersist pos(ntpolyname);
417 cout << " SimpleOffsetEstimator::run()/Info : Writing poly ntuple to PPF file " << ntpolyname << endl;
418 pos << xntp;
419 }
420}
Note: See TracBrowser for help on using the repository browser.