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

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

Conversion XNTuple en NTuple ds simoffset.cc - Reza 31/5/2002

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