1 | //
|
---|
2 | // $Id: histos.cc,v 1.1.1.1 1999-04-09 17:57:56 ansari Exp $
|
---|
3 | //
|
---|
4 |
|
---|
5 | #include "defs.h"
|
---|
6 | #include <string.h>
|
---|
7 | #include <stdio.h>
|
---|
8 | #include <math.h>
|
---|
9 | #include "histos.h"
|
---|
10 | #include "perrors.h"
|
---|
11 | #include "cvector.h"
|
---|
12 | #include "poly.h"
|
---|
13 | #include "strutil.h"
|
---|
14 | #include "generalfit.h"
|
---|
15 |
|
---|
16 | //++
|
---|
17 | // Class Histo
|
---|
18 | // Lib Outils++
|
---|
19 | // include histos.h
|
---|
20 | //
|
---|
21 | // Classe d'histogrammes 1D
|
---|
22 | //--
|
---|
23 |
|
---|
24 | //++
|
---|
25 | // Titre Constructeurs
|
---|
26 | //--
|
---|
27 |
|
---|
28 | /********* Methode *********/
|
---|
29 | //++
|
---|
30 | Histo::Histo()
|
---|
31 | //
|
---|
32 | //--
|
---|
33 | : data(NULL), err2(NULL),
|
---|
34 | under(0), over(0), nHist(0), nEntries(0),
|
---|
35 | bins(0), min(0), max(0),
|
---|
36 | binWidth(0)
|
---|
37 | {
|
---|
38 | END_CONSTRUCTOR
|
---|
39 | }
|
---|
40 |
|
---|
41 | /********* Methode *********/
|
---|
42 | //++
|
---|
43 | Histo::Histo(float xMin, float xMax, int nBin)
|
---|
44 | //
|
---|
45 | //--
|
---|
46 | : data(new float[nBin]), err2(NULL),
|
---|
47 | under(0), over(0), nHist(0), nEntries(0),
|
---|
48 | bins(nBin), min(xMin), max(xMax),
|
---|
49 | binWidth((max - min)/nBin)
|
---|
50 | {
|
---|
51 | Zero();
|
---|
52 | END_CONSTRUCTOR
|
---|
53 | }
|
---|
54 |
|
---|
55 | /********* Methode *********/
|
---|
56 | //++
|
---|
57 | Histo::Histo(char *flnm)
|
---|
58 | //
|
---|
59 | //--
|
---|
60 | : data(NULL), err2(NULL),
|
---|
61 | under(0), over(0), nHist(0), nEntries(0),
|
---|
62 | bins(0), min(0), max(0), binWidth(0)
|
---|
63 | {
|
---|
64 | PInPersist s(flnm);
|
---|
65 | Read(s);
|
---|
66 | END_CONSTRUCTOR
|
---|
67 | }
|
---|
68 |
|
---|
69 | /********* Methode *********/
|
---|
70 | //++
|
---|
71 | Histo::Histo(const Histo& H)
|
---|
72 | //
|
---|
73 | //--
|
---|
74 | : data(new float[H.bins]), err2(NULL),
|
---|
75 | under(H.under), over(H.over), nHist(H.nHist), nEntries(H.nEntries),
|
---|
76 | bins(H.bins), min(H.min), max(H.max),
|
---|
77 | binWidth(H.binWidth)
|
---|
78 | {
|
---|
79 | memcpy(data, H.data, bins*sizeof(float));
|
---|
80 | if( H.err2 != NULL ) {
|
---|
81 | err2 = new double[bins];
|
---|
82 | memcpy(err2, H.err2, bins*sizeof(double));
|
---|
83 | }
|
---|
84 | END_CONSTRUCTOR
|
---|
85 | }
|
---|
86 |
|
---|
87 | /********* Methode *********/
|
---|
88 | void Histo::Delete()
|
---|
89 | {
|
---|
90 | if( data != NULL ) { delete[] data; data = NULL;}
|
---|
91 | if( err2 != NULL ) { delete[] err2; err2 = NULL;}
|
---|
92 | under = over = min = max = binWidth= 0.;
|
---|
93 | nHist = 0.;
|
---|
94 | nEntries = bins = 0;
|
---|
95 | }
|
---|
96 |
|
---|
97 | /********* Methode *********/
|
---|
98 | Histo::~Histo()
|
---|
99 | {
|
---|
100 | Delete();
|
---|
101 | }
|
---|
102 |
|
---|
103 | //++
|
---|
104 | // Titre Methodes
|
---|
105 | //--
|
---|
106 |
|
---|
107 | /********* Methode *********/
|
---|
108 | //++
|
---|
109 | void Histo::Zero()
|
---|
110 | //
|
---|
111 | // Remise a zero
|
---|
112 | //--
|
---|
113 | {
|
---|
114 | memset(data, 0, bins*sizeof(float));
|
---|
115 | under = over = 0;
|
---|
116 | nHist = 0;
|
---|
117 | nEntries = 0;
|
---|
118 | if( err2 != NULL ) memset(err2, 0, bins*sizeof(double));
|
---|
119 | }
|
---|
120 |
|
---|
121 | /********* Methode *********/
|
---|
122 | //++
|
---|
123 | void Histo::Errors()
|
---|
124 | //
|
---|
125 | // Pour avoir le calcul des erreurs
|
---|
126 | //--
|
---|
127 | {
|
---|
128 | if( bins > 0 ) {
|
---|
129 | if(err2==NULL) err2 = new double[bins];
|
---|
130 | memset(err2, 0, bins*sizeof(double));
|
---|
131 | }
|
---|
132 | }
|
---|
133 |
|
---|
134 | /********* Methode *********/
|
---|
135 | //++
|
---|
136 | Histo& Histo::operator = (const Histo& h)
|
---|
137 | //
|
---|
138 | //--
|
---|
139 | {
|
---|
140 | if(this == &h) return *this;
|
---|
141 | if( h.bins > bins ) Delete();
|
---|
142 | if(!data) data = new float[h.bins];
|
---|
143 | if( !h.err2 && err2 ) { delete [] err2; err2=NULL;}
|
---|
144 | if( h.err2 && !err2 ) err2 = new double[h.bins];
|
---|
145 |
|
---|
146 | under = h.under;
|
---|
147 | over = h.over;
|
---|
148 | nHist = h.nHist;
|
---|
149 | nEntries = h.nEntries;
|
---|
150 | bins = h.bins;
|
---|
151 | min = h.min;
|
---|
152 | max = h.max;
|
---|
153 | binWidth = h.binWidth;
|
---|
154 |
|
---|
155 | memcpy(data, h.data, bins*sizeof(float));
|
---|
156 | if(err2) memcpy(err2, h.err2, bins*sizeof(double));
|
---|
157 |
|
---|
158 | return *this;
|
---|
159 | }
|
---|
160 |
|
---|
161 | /********* Methode *********/
|
---|
162 | //++
|
---|
163 | Histo& Histo::operator *= (double b)
|
---|
164 | //
|
---|
165 | //--
|
---|
166 | {
|
---|
167 | double b2 = b*b;
|
---|
168 | for(int i=0;i<bins;i++) {
|
---|
169 | data[i] *= b;
|
---|
170 | if(err2) err2[i] *= b2;
|
---|
171 | }
|
---|
172 | under *= b;
|
---|
173 | over *= b;
|
---|
174 | nHist *= b;
|
---|
175 |
|
---|
176 | return *this;
|
---|
177 | }
|
---|
178 |
|
---|
179 | //++
|
---|
180 | Histo& Histo::operator /= (double b)
|
---|
181 | //
|
---|
182 | //--
|
---|
183 | {
|
---|
184 | if (b==0.) THROW(inconsistentErr);
|
---|
185 | double b2 = b*b;
|
---|
186 | for(int i=0;i<bins;i++) {
|
---|
187 | data[i] /= b;
|
---|
188 | if(err2) err2[i] /= b2;
|
---|
189 | }
|
---|
190 | under /= b;
|
---|
191 | over /= b;
|
---|
192 | nHist /= b;
|
---|
193 |
|
---|
194 | return *this;
|
---|
195 | }
|
---|
196 |
|
---|
197 | //++
|
---|
198 | Histo& Histo::operator += (double b)
|
---|
199 | //
|
---|
200 | //--
|
---|
201 | {
|
---|
202 | for(int i=0;i<bins;i++) data[i] += b;
|
---|
203 | under += b;
|
---|
204 | over += b;
|
---|
205 | nHist += bins*b;
|
---|
206 |
|
---|
207 | return *this;
|
---|
208 | }
|
---|
209 |
|
---|
210 | //++
|
---|
211 | Histo& Histo::operator -= (double b)
|
---|
212 | //
|
---|
213 | //--
|
---|
214 | {
|
---|
215 | for(int i=0;i<bins;i++) data[i] -= b;
|
---|
216 | under -= b;
|
---|
217 | over -= b;
|
---|
218 | nHist -= bins*b;
|
---|
219 |
|
---|
220 | return *this;
|
---|
221 | }
|
---|
222 |
|
---|
223 | /********* Methode *********/
|
---|
224 | //++
|
---|
225 | Histo operator * (const Histo& a, double b)
|
---|
226 | //
|
---|
227 | //--
|
---|
228 | {
|
---|
229 | Histo result(a);
|
---|
230 | return (result *= b);
|
---|
231 | }
|
---|
232 |
|
---|
233 | //++
|
---|
234 | Histo operator * (double b, const Histo& a)
|
---|
235 | //
|
---|
236 | //--
|
---|
237 | {
|
---|
238 | Histo result(a);
|
---|
239 | return (result *= b);
|
---|
240 | }
|
---|
241 |
|
---|
242 | //++
|
---|
243 | Histo operator / (const Histo& a, double b)
|
---|
244 | //
|
---|
245 | //--
|
---|
246 | {
|
---|
247 | Histo result(a);
|
---|
248 | return (result /= b);
|
---|
249 | }
|
---|
250 |
|
---|
251 | //++
|
---|
252 | Histo operator + (const Histo& a, double b)
|
---|
253 | //
|
---|
254 | //--
|
---|
255 | {
|
---|
256 | Histo result(a);
|
---|
257 | return (result += b);
|
---|
258 | }
|
---|
259 |
|
---|
260 | //++
|
---|
261 | Histo operator + (double b, const Histo& a)
|
---|
262 | //
|
---|
263 | //--
|
---|
264 | {
|
---|
265 | Histo result(a);
|
---|
266 | return (result += b);
|
---|
267 | }
|
---|
268 |
|
---|
269 | //++
|
---|
270 | Histo operator - (const Histo& a, double b)
|
---|
271 | //
|
---|
272 | //--
|
---|
273 | {
|
---|
274 | Histo result(a);
|
---|
275 | return (result -= b);
|
---|
276 | }
|
---|
277 |
|
---|
278 | //++
|
---|
279 | Histo operator - (double b, const Histo& a)
|
---|
280 | //
|
---|
281 | //--
|
---|
282 | {
|
---|
283 | Histo result(a);
|
---|
284 | result *= -1.;
|
---|
285 | return (result += b);
|
---|
286 | }
|
---|
287 |
|
---|
288 | /********* Methode *********/
|
---|
289 | //++
|
---|
290 | Histo& Histo::operator += (const Histo& a)
|
---|
291 | //
|
---|
292 | //--
|
---|
293 | {
|
---|
294 | if(bins!=a.bins) THROW(sizeMismatchErr);
|
---|
295 | for(int i=0;i<bins;i++) {
|
---|
296 | data[i] += a(i);
|
---|
297 | if(err2 && a.err2) err2[i] += a.Error2(i);
|
---|
298 | else err2[i] = 0.;
|
---|
299 | }
|
---|
300 | under += a.under;
|
---|
301 | over += a.over;
|
---|
302 | nHist += a.nHist;
|
---|
303 | nEntries += a.nEntries;
|
---|
304 |
|
---|
305 | return *this;
|
---|
306 | }
|
---|
307 |
|
---|
308 | //++
|
---|
309 | Histo& Histo::operator -= (const Histo& a)
|
---|
310 | //
|
---|
311 | //--
|
---|
312 | {
|
---|
313 | if(bins!=a.bins) THROW(sizeMismatchErr);
|
---|
314 | for(int i=0;i<bins;i++) {
|
---|
315 | data[i] -= a(i);
|
---|
316 | if(err2 && a.err2) err2[i] += a.Error2(i);
|
---|
317 | else err2[i] = 0.;
|
---|
318 | }
|
---|
319 | under -= a.under;
|
---|
320 | over -= a.over;
|
---|
321 | nHist -= a.nHist;
|
---|
322 | nEntries += a.nEntries;
|
---|
323 |
|
---|
324 | return *this;
|
---|
325 | }
|
---|
326 |
|
---|
327 | //++
|
---|
328 | Histo& Histo::operator *= (const Histo& a)
|
---|
329 | //
|
---|
330 | //--
|
---|
331 | {
|
---|
332 | if(bins!=a.bins) THROW(sizeMismatchErr);
|
---|
333 | nHist = 0.;
|
---|
334 | for(int i=0;i<bins;i++) {
|
---|
335 | if(err2 && a.err2)
|
---|
336 | err2[i] = a(i)*a(i)*err2[i] + data[i]*data[i]*a.Error2(i);
|
---|
337 | else err2[i] = 0.;
|
---|
338 | data[i] *= a(i);
|
---|
339 | nHist += data[i];
|
---|
340 | }
|
---|
341 | under *= a.under;
|
---|
342 | over *= a.over;
|
---|
343 | nEntries += a.nEntries;
|
---|
344 |
|
---|
345 | return *this;
|
---|
346 | }
|
---|
347 |
|
---|
348 | //++
|
---|
349 | Histo& Histo::operator /= (const Histo& a)
|
---|
350 | //
|
---|
351 | //--
|
---|
352 | {
|
---|
353 | if(bins!=a.bins) THROW(sizeMismatchErr);
|
---|
354 | nHist = 0.;
|
---|
355 | for(int i=0;i<bins;i++) {
|
---|
356 | if(a(i)==0.) {
|
---|
357 | data[i]=0.;
|
---|
358 | if(err2) err2[i]=0.;
|
---|
359 | continue;
|
---|
360 | }
|
---|
361 | if(err2 && a.err2)
|
---|
362 | err2[i] = (err2[i] + data[i]/a(i)*data[i]/a(i)*a.Error2(i))
|
---|
363 | /(a(i)*a(i));
|
---|
364 | else err2[i] = 0.;
|
---|
365 | data[i] /= a(i);
|
---|
366 | nHist += data[i];
|
---|
367 | }
|
---|
368 | if(a.under!=0.) under /= a.under; else under = 0.;
|
---|
369 | if(a.over!=0.) over /= a.over; else over = 0.;
|
---|
370 | nEntries += a.nEntries;
|
---|
371 |
|
---|
372 | return *this;
|
---|
373 | }
|
---|
374 |
|
---|
375 | /********* Methode *********/
|
---|
376 | //++
|
---|
377 | Histo operator + (const Histo& a, const Histo& b)
|
---|
378 | //
|
---|
379 | //--
|
---|
380 | {
|
---|
381 | if (b.bins!=a.bins) THROW(sizeMismatchErr);
|
---|
382 | Histo c(a);
|
---|
383 | return (c += b);
|
---|
384 | }
|
---|
385 |
|
---|
386 | //++
|
---|
387 | Histo operator - (const Histo& a, const Histo& b)
|
---|
388 | //
|
---|
389 | //--
|
---|
390 | {
|
---|
391 | if (b.bins!=a.bins) THROW(sizeMismatchErr);
|
---|
392 | Histo c(a);
|
---|
393 | return (c -= b);
|
---|
394 | }
|
---|
395 |
|
---|
396 | //++
|
---|
397 | Histo operator * (const Histo& a, const Histo& b)
|
---|
398 | //
|
---|
399 | //--
|
---|
400 | {
|
---|
401 | if (b.bins!=a.bins) THROW(sizeMismatchErr);
|
---|
402 | Histo c(a);
|
---|
403 | return (c *= b);
|
---|
404 | }
|
---|
405 |
|
---|
406 | //++
|
---|
407 | Histo operator / (const Histo& a, const Histo& b)
|
---|
408 | //
|
---|
409 | //--
|
---|
410 | {
|
---|
411 | if (b.bins!=a.bins) THROW(sizeMismatchErr);
|
---|
412 | Histo c(a);
|
---|
413 | return (c /= b);
|
---|
414 | }
|
---|
415 |
|
---|
416 | /********* Methode *********/
|
---|
417 | //++
|
---|
418 | void Histo::GetAbsc(Vector &v)
|
---|
419 | //
|
---|
420 | // Remplissage d'un tableau avec la valeur des abscisses
|
---|
421 | //--
|
---|
422 | {
|
---|
423 | v.Realloc(bins);
|
---|
424 | for(int i=0;i<bins;i++) v(i) = BinLowEdge(i);
|
---|
425 | return;
|
---|
426 | }
|
---|
427 |
|
---|
428 | //++
|
---|
429 | void Histo::GetValue(Vector &v)
|
---|
430 | //
|
---|
431 | // Remplissage d'un tableau avec la valeur du contenu
|
---|
432 | //--
|
---|
433 | {
|
---|
434 | v.Realloc(bins);
|
---|
435 | for(int i=0;i<bins;i++) v(i) = data[i];
|
---|
436 | return;
|
---|
437 | }
|
---|
438 |
|
---|
439 | //++
|
---|
440 | void Histo::GetError2(Vector &v)
|
---|
441 | //
|
---|
442 | // Remplissage d'un tableau avec la valeur des erreurs au carre
|
---|
443 | //--
|
---|
444 | {
|
---|
445 | v.Realloc(bins);
|
---|
446 | if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
|
---|
447 | for(int i=0;i<bins;i++) v(i) = err2[i];
|
---|
448 | return;
|
---|
449 | }
|
---|
450 |
|
---|
451 | //++
|
---|
452 | void Histo::GetError(Vector &v)
|
---|
453 | //
|
---|
454 | // Remplissage d'un tableau avec la valeur des erreurs
|
---|
455 | //--
|
---|
456 | {
|
---|
457 | v.Realloc(bins);
|
---|
458 | if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
|
---|
459 | for(int i=0;i<bins;i++) v(i) = Error(i);
|
---|
460 | return;
|
---|
461 | }
|
---|
462 |
|
---|
463 | /********* Methode *********/
|
---|
464 | //++
|
---|
465 | void Histo::PutValue(Vector &v, int ierr)
|
---|
466 | //
|
---|
467 | // Remplissage du contenu de l'histo avec les valeurs d'un vecteur
|
---|
468 | //--
|
---|
469 | {
|
---|
470 | if(v.NElts()<bins) THROW(sizeMismatchErr);
|
---|
471 | for(int i=0;i<bins;i++) {
|
---|
472 | data[i] = v(i);
|
---|
473 | if(err2&&ierr) err2[i] = fabs(v(i));
|
---|
474 | }
|
---|
475 | return;
|
---|
476 | }
|
---|
477 |
|
---|
478 | //++
|
---|
479 | void Histo::PutValueAdd(Vector &v, int ierr)
|
---|
480 | //
|
---|
481 | // Addition du contenu de l'histo avec les valeurs d'un vecteur
|
---|
482 | //--
|
---|
483 | {
|
---|
484 | if(v.NElts()<bins) THROW(sizeMismatchErr);
|
---|
485 | for(int i=0;i<bins;i++) {
|
---|
486 | data[i] += v(i);
|
---|
487 | if(err2&&ierr) err2[i] += fabs(v(i));
|
---|
488 | }
|
---|
489 | return;
|
---|
490 | }
|
---|
491 |
|
---|
492 | //++
|
---|
493 | void Histo::PutError2(Vector &v)
|
---|
494 | //
|
---|
495 | // Remplissage des erreurs au carre de l'histo avec les valeurs d'un vecteur
|
---|
496 | //--
|
---|
497 | {
|
---|
498 | if(v.NElts()<bins) THROW(sizeMismatchErr);
|
---|
499 | if(!err2) Errors();
|
---|
500 | for(int i=0;i<bins;i++) err2[i] = v(i);
|
---|
501 | return;
|
---|
502 | }
|
---|
503 |
|
---|
504 | //++
|
---|
505 | void Histo::PutError2Add(Vector &v)
|
---|
506 | //
|
---|
507 | // Addition des erreurs au carre de l'histo avec les valeurs d'un vecteur
|
---|
508 | //--
|
---|
509 | {
|
---|
510 | if(v.NElts()<bins) THROW(sizeMismatchErr);
|
---|
511 | if(!err2) Errors();
|
---|
512 | for(int i=0;i<bins;i++) if(v(i)>0.) err2[i] += v(i);
|
---|
513 | return;
|
---|
514 | }
|
---|
515 |
|
---|
516 | //++
|
---|
517 | void Histo::PutError(Vector &v)
|
---|
518 | //
|
---|
519 | // Remplissage des erreurs de l'histo avec les valeurs d'un vecteur
|
---|
520 | //--
|
---|
521 | {
|
---|
522 | if(v.NElts()<bins) THROW(sizeMismatchErr);
|
---|
523 | if(!err2) Errors();
|
---|
524 | for(int i=0;i<bins;i++)
|
---|
525 | if(v(i)>0.) err2[i]=v(i)*v(i); else err2[i]=-v(i)*v(i);
|
---|
526 | return;
|
---|
527 | }
|
---|
528 |
|
---|
529 | /********* Methode *********/
|
---|
530 | //++
|
---|
531 | void Histo::Add(float x, float w)
|
---|
532 | //
|
---|
533 | // Addition du contenu de l'histo pour abscisse x poids w
|
---|
534 | //--
|
---|
535 | {
|
---|
536 | int numBin = FindBin(x);
|
---|
537 | if (numBin<0) under += w;
|
---|
538 | else if (numBin>=bins) over += w;
|
---|
539 | else {
|
---|
540 | data[numBin] += w;
|
---|
541 | if(err2!=NULL) err2[numBin] += w*w;
|
---|
542 | nHist += w;
|
---|
543 | nEntries++;
|
---|
544 | }
|
---|
545 | }
|
---|
546 |
|
---|
547 | /********* Methode *********/
|
---|
548 | //++
|
---|
549 | void Histo::AddBin(int numBin, float w)
|
---|
550 | //
|
---|
551 | // Addition du contenu de l'histo bin numBin poids w
|
---|
552 | //--
|
---|
553 | {
|
---|
554 | if (numBin<0) under += w;
|
---|
555 | else if (numBin>=bins) over += w;
|
---|
556 | else {
|
---|
557 | data[numBin] += w;
|
---|
558 | if(err2!=NULL) err2[numBin] += w*w;
|
---|
559 | nHist += w;
|
---|
560 | nEntries++;
|
---|
561 | }
|
---|
562 | }
|
---|
563 |
|
---|
564 | /********* Methode *********/
|
---|
565 | //++
|
---|
566 | void Histo::SetBin(float x, float w)
|
---|
567 | //
|
---|
568 | // Remplissage du contenu de l'histo pour abscisse x poids w
|
---|
569 | //--
|
---|
570 | {
|
---|
571 | int numBin = FindBin(x);
|
---|
572 | SetBin(numBin,w);
|
---|
573 | }
|
---|
574 |
|
---|
575 | /********* Methode *********/
|
---|
576 | //++
|
---|
577 | void Histo::SetBin(int numBin, float w)
|
---|
578 | //
|
---|
579 | // Remplissage du contenu de l'histo pour numBin poids w
|
---|
580 | //--
|
---|
581 | {
|
---|
582 | if (numBin<0) under = w;
|
---|
583 | else if (numBin>=bins) over = w;
|
---|
584 | else {
|
---|
585 | nHist -= data[numBin];
|
---|
586 | data[numBin] = w;
|
---|
587 | nHist += w;
|
---|
588 | }
|
---|
589 | }
|
---|
590 |
|
---|
591 | /********* Methode *********/
|
---|
592 | //++
|
---|
593 | void Histo::SetErr2(float x, double e2)
|
---|
594 | //
|
---|
595 | // Remplissage des erreurs au carre pour abscisse x
|
---|
596 | //--
|
---|
597 | {
|
---|
598 | int numBin = FindBin(x);
|
---|
599 | SetErr2(numBin,e2);
|
---|
600 | }
|
---|
601 |
|
---|
602 | /********* Methode *********/
|
---|
603 | //++
|
---|
604 | void Histo::SetErr2(int numBin, double e2)
|
---|
605 | //
|
---|
606 | // Remplissage des erreurs au carre pour numBin poids
|
---|
607 | //--
|
---|
608 | {
|
---|
609 | if( err2==NULL) return;
|
---|
610 | if ( numBin<0 || numBin>=bins ) return;
|
---|
611 | err2[numBin] = e2;
|
---|
612 | }
|
---|
613 |
|
---|
614 | /********* Methode *********/
|
---|
615 | //++
|
---|
616 | void Histo::SetErr(float x, float e)
|
---|
617 | //
|
---|
618 | // Remplissage des erreurs pour abscisse x
|
---|
619 | //--
|
---|
620 | {
|
---|
621 | SetErr2(x, (double) e*e);
|
---|
622 | }
|
---|
623 |
|
---|
624 | /********* Methode *********/
|
---|
625 | //++
|
---|
626 | void Histo::SetErr(int numBin, float e)
|
---|
627 | //
|
---|
628 | // Remplissage des erreurs pour numBin
|
---|
629 | //--
|
---|
630 | {
|
---|
631 | SetErr2(numBin, (double) e*e);
|
---|
632 | }
|
---|
633 |
|
---|
634 | /********* Methode *********/
|
---|
635 | //++
|
---|
636 | int Histo::IMax() const
|
---|
637 | //
|
---|
638 | // Numero du bin ayant le contenu maximum
|
---|
639 | //--
|
---|
640 | {
|
---|
641 | int imx=0;
|
---|
642 | if(bins==1) return imx;
|
---|
643 | float mx=data[0];
|
---|
644 | for (int i=1; i<bins; i++)
|
---|
645 | if (data[i]>mx) {imx = i; mx=data[i];}
|
---|
646 | return imx;
|
---|
647 | }
|
---|
648 |
|
---|
649 | /********* Methode *********/
|
---|
650 | //++
|
---|
651 | int Histo::IMin() const
|
---|
652 | //
|
---|
653 | // Numero du bin ayant le contenu minimum
|
---|
654 | //--
|
---|
655 | {
|
---|
656 | int imx=0;
|
---|
657 | if(bins==1) return imx;
|
---|
658 | float mx=data[0];
|
---|
659 | for (int i=1; i<bins; i++)
|
---|
660 | if (data[i]<mx) {imx = i; mx=data[i];}
|
---|
661 | return imx;
|
---|
662 | }
|
---|
663 |
|
---|
664 | /********* Methode *********/
|
---|
665 | //++
|
---|
666 | float Histo::VMax() const
|
---|
667 | //
|
---|
668 | // Valeur le contenu maximum
|
---|
669 | //--
|
---|
670 | {
|
---|
671 | float mx=data[0];
|
---|
672 | if(bins==1) return mx;
|
---|
673 | for (int i=1;i<bins;i++) if(data[i]>mx) mx=data[i];
|
---|
674 | return mx;
|
---|
675 | }
|
---|
676 |
|
---|
677 | /********* Methode *********/
|
---|
678 | //++
|
---|
679 | float Histo::VMin() const
|
---|
680 | //
|
---|
681 | // Valeur le contenu minimum
|
---|
682 | //--
|
---|
683 | {
|
---|
684 | float mx=data[0];
|
---|
685 | if(bins==1) return mx;
|
---|
686 | for (int i=1;i<bins;i++) if(data[i]<mx) mx=data[i];
|
---|
687 | return mx;
|
---|
688 | }
|
---|
689 |
|
---|
690 | /********* Methode *********/
|
---|
691 | //++
|
---|
692 | float Histo::Mean() const
|
---|
693 | //
|
---|
694 | // Valeur moyenne
|
---|
695 | //--
|
---|
696 | {
|
---|
697 | double n = 0;
|
---|
698 | double sx = 0;
|
---|
699 | for (int i=0; i<bins; i++) {
|
---|
700 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
701 | n += v;
|
---|
702 | sx += BinCenter(i)*v;
|
---|
703 | }
|
---|
704 | if(n>0.) return sx/n; else return min;
|
---|
705 | }
|
---|
706 |
|
---|
707 | /********* Methode *********/
|
---|
708 | //++
|
---|
709 | float Histo::Sigma() const
|
---|
710 | //
|
---|
711 | // Valeur du sigma
|
---|
712 | //--
|
---|
713 | {
|
---|
714 | double n = 0;
|
---|
715 | double sx = 0;
|
---|
716 | double sx2= 0;
|
---|
717 | for (int i=0; i<bins; i++) {
|
---|
718 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
719 | n += v;
|
---|
720 | sx += BinCenter(i)*v;
|
---|
721 | sx2+= BinCenter(i)*BinCenter(i)*v;
|
---|
722 | }
|
---|
723 | // attention a l'erreur d'arrondi si un seul bin rempli!!
|
---|
724 | if( n == 0 ) return 0.;
|
---|
725 | sx2 = sx2/n - (sx/n)*(sx/n);
|
---|
726 | if( sx2>0. ) return sqrt( sx2 );
|
---|
727 | else return 0.;
|
---|
728 | }
|
---|
729 |
|
---|
730 | /********* Methode *********/
|
---|
731 | //++
|
---|
732 | float Histo::MeanLH(int il,int ih) const
|
---|
733 | //
|
---|
734 | // Valeur de la moyenne calculee entre les bin il et ih
|
---|
735 | //--
|
---|
736 | {
|
---|
737 | if( ih < il ) { il = 0; ih = bins-1; }
|
---|
738 | if( il < 0 ) il = 0;
|
---|
739 | if( ih >= bins ) ih = bins-1;
|
---|
740 | double n = 0;
|
---|
741 | double sx = 0;
|
---|
742 | for (int i=il; i<=ih; i++) {
|
---|
743 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
744 | n += v;
|
---|
745 | sx += BinCenter(i)*v;
|
---|
746 | }
|
---|
747 | if(n>0.) return sx/n; else return BinLowEdge(il);
|
---|
748 | }
|
---|
749 |
|
---|
750 | /********* Methode *********/
|
---|
751 | //++
|
---|
752 | float Histo::SigmaLH(int il,int ih) const
|
---|
753 | //
|
---|
754 | // Valeur de la moyenne calculee entre les bin il et ih
|
---|
755 | //--
|
---|
756 | {
|
---|
757 | if( ih < il ) { il = 0; ih = bins - 1; }
|
---|
758 | if( il < 0 ) il = 0;
|
---|
759 | if( ih >= bins ) ih = bins - 1;
|
---|
760 | double n = 0;
|
---|
761 | double sx = 0;
|
---|
762 | double sx2= 0;
|
---|
763 | for (int i=il; i<=ih; i++) {
|
---|
764 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
765 | n += v;
|
---|
766 | sx += BinCenter(i)*v;
|
---|
767 | sx2+= BinCenter(i)*BinCenter(i)*v;
|
---|
768 | }
|
---|
769 | if( n == 0 ) return 0.;
|
---|
770 | sx2 = sx2/n - (sx/n)*(sx/n);
|
---|
771 | if( sx2>0. ) return sqrt( sx2 );
|
---|
772 | else return 0.;
|
---|
773 | }
|
---|
774 |
|
---|
775 | /********* Methode *********/
|
---|
776 | //++
|
---|
777 | float Histo::Mean(float x0, float dx) const
|
---|
778 | //
|
---|
779 | // Valeur de la moyenne calculee entre x0-dx et x0+dx
|
---|
780 | //--
|
---|
781 | {
|
---|
782 | double sdata = 0;
|
---|
783 | double sx = 0;
|
---|
784 | int iMin = FindBin(x0-dx);
|
---|
785 | if (iMin<0) iMin = 0;
|
---|
786 | int iMax = FindBin(x0+dx);
|
---|
787 | if (iMax>bins) iMax=bins;
|
---|
788 | for (int i=iMin; i<iMax; i++) {
|
---|
789 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
790 | sx += BinCenter(i)*v;
|
---|
791 | sdata += v;
|
---|
792 | }
|
---|
793 | if(sdata>0.) return sx/sdata; else return min;
|
---|
794 | }
|
---|
795 |
|
---|
796 | /********* Methode *********/
|
---|
797 | //++
|
---|
798 | float Histo::Sigma(float x0, float dx) const
|
---|
799 | //
|
---|
800 | // Valeur du sigma calcule entre x0-dx et x0+dx
|
---|
801 | //--
|
---|
802 | {
|
---|
803 | double sx = 0;
|
---|
804 | double sx2= 0;
|
---|
805 | double sdata = 0;
|
---|
806 | int iMin = FindBin(x0-dx);
|
---|
807 | if (iMin<0) iMin = 0;
|
---|
808 | int iMax = FindBin(x0+dx);
|
---|
809 | if (iMax>bins) iMax=bins;
|
---|
810 | for (int i=iMin; i<iMax; i++) {
|
---|
811 | double v = (data[i]>=0.) ? data[i] : -data[i];
|
---|
812 | sx += BinCenter(i)*v;
|
---|
813 | sx2+= BinCenter(i)*BinCenter(i)*v;
|
---|
814 | sdata += v;
|
---|
815 | }
|
---|
816 | if(sdata>0.) return sqrt( sx2/sdata - (sx/sdata)*(sx/sdata) );
|
---|
817 | else return 0.;
|
---|
818 | }
|
---|
819 |
|
---|
820 | /********* Methode *********/
|
---|
821 | //++
|
---|
822 | float Histo::CleanedMean(float& sigma) const
|
---|
823 | //
|
---|
824 | // Valeur de la moyenne et du sigma nettoyes
|
---|
825 | //--
|
---|
826 | {
|
---|
827 | if (!nHist) return 0;
|
---|
828 | // on fait ca de facon bete, on pourra raffiner apres...
|
---|
829 | float x0 = Mean();
|
---|
830 | float s = Sigma()+binWidth;
|
---|
831 |
|
---|
832 | for (int i=0; i<3; i++) {
|
---|
833 | float xx0 = Mean(x0,5*s);
|
---|
834 | s = Sigma(x0,5*s)+binWidth;
|
---|
835 | x0 = xx0;
|
---|
836 | }
|
---|
837 | sigma = s;
|
---|
838 | return x0;
|
---|
839 | }
|
---|
840 |
|
---|
841 | /********* Methode *********/
|
---|
842 | //++
|
---|
843 | float Histo::CleanedMean() const
|
---|
844 | //
|
---|
845 | // Valeur de la moyenne nettoyee
|
---|
846 | //--
|
---|
847 | {
|
---|
848 | if (!nHist) return 0;
|
---|
849 | float s=0;
|
---|
850 | return CleanedMean(s);
|
---|
851 | }
|
---|
852 |
|
---|
853 | /********* Methode *********/
|
---|
854 | //++
|
---|
855 | int Histo::BinNonNul() const
|
---|
856 | //
|
---|
857 | // Retourne le nombre de bins non-nul
|
---|
858 | //--
|
---|
859 | {
|
---|
860 | int non=0;
|
---|
861 | for (int i=0;i<bins;i++) if( data[i] != 0. ) non++;
|
---|
862 | return non;
|
---|
863 | }
|
---|
864 |
|
---|
865 | /********* Methode *********/
|
---|
866 | //++
|
---|
867 | int Histo::ErrNonNul() const
|
---|
868 | //
|
---|
869 | // Retourne le nombre de bins ayant une erreur non-nulle
|
---|
870 | //--
|
---|
871 | {
|
---|
872 | if(err2==NULL) return -1;
|
---|
873 | int non=0;
|
---|
874 | for (int i=0;i<bins;i++) if( err2[i] != 0. ) non++;
|
---|
875 | return non;
|
---|
876 | }
|
---|
877 |
|
---|
878 | /********* Methode *********/
|
---|
879 | //++
|
---|
880 | int Histo::BinPercent(float per) const
|
---|
881 | //
|
---|
882 | // Renvoie le numero de bin tel que il y ait "per" pourcent entrees
|
---|
883 | // entre le bin 0 et ce bin (y compris le contenu de ce bin).
|
---|
884 | //--
|
---|
885 | {
|
---|
886 | double n = nHist*per;
|
---|
887 | double s = 0.;
|
---|
888 | int i;
|
---|
889 | for(i=0; i<bins; i++ ) {
|
---|
890 | s += data[i];
|
---|
891 | if( s >= n ) break;
|
---|
892 | }
|
---|
893 | if( i == bins ) i = bins-1;
|
---|
894 | return i;
|
---|
895 | }
|
---|
896 |
|
---|
897 | /********* Methode *********/
|
---|
898 | //++
|
---|
899 | int Histo::BinPercent(float x,float per,int& imin,int& imax)
|
---|
900 | //
|
---|
901 | // Renvoie les numeros de bins imin,imax (0=<i<bins) tels que:
|
---|
902 | //| notons I = bin contenant l'abscisse x
|
---|
903 | //| N1 = nombre d'entrees entre bin 0 et I compris
|
---|
904 | //| N2 = nombre d'entrees entre bin I et bins-1 compris
|
---|
905 | //| imin = bin tel que nombre d'entrees entre imin et I = N1 * per
|
---|
906 | //| imax = bin tel que nombre d'entrees entre I et imax = N2 * per
|
---|
907 | //| Return: <0 si echec
|
---|
908 | //| min(I-imin,imax-I) si ok
|
---|
909 | //--
|
---|
910 | {
|
---|
911 | imin = imax = -1;
|
---|
912 | if( per <= 0. ) return -1;
|
---|
913 |
|
---|
914 | int I = FindBin(x);
|
---|
915 | if( I<0 || I>=bins ) return -2;
|
---|
916 |
|
---|
917 | double N1 = 0.; for(int i=0; i<=I; i++) N1 += data[i]; N1 *= per;
|
---|
918 | double N2 = 0.; {for(int i=I; i<bins; i++) N2 += data[i]; N2 *= per;}
|
---|
919 |
|
---|
920 | double n = 0.;
|
---|
921 | for(imin=I; imin>=0; imin--) { n += data[imin]; if(n>=N1) break; }
|
---|
922 | if( imin<0 ) imin = 0;
|
---|
923 | // cout<<"I="<<I<<" imin="<<imin<<" n="<<n<<" N1="<<N1<<endl;
|
---|
924 |
|
---|
925 | n = 0.;
|
---|
926 | for(imax=I; imax<bins; imax++) { n += data[imax]; if(n>=N2) break; }
|
---|
927 | if( imax>=bins ) imax = bins-1;
|
---|
928 | // cout<<"I="<<I<<" imax="<<imax<<" n="<<n<<" N2="<<N2<<endl;
|
---|
929 |
|
---|
930 | return ( imax-I > I-imin ) ? I-imin : imax-I ;
|
---|
931 | }
|
---|
932 |
|
---|
933 | /********* Methode *********/
|
---|
934 | //++
|
---|
935 | int Histo::BinPercent(float x,float per,float& xmin,float& xmax)
|
---|
936 | //
|
---|
937 | // Idem precedent mais renvoie xmin et xmax
|
---|
938 | //--
|
---|
939 | {
|
---|
940 | xmin = xmax = 0.;
|
---|
941 | int imin,imax;
|
---|
942 | int rc = BinPercent(x,per,imin,imax);
|
---|
943 | if( rc >= 0 ) { xmin = BinLowEdge(imin); xmax = BinHighEdge(imax); }
|
---|
944 | return rc;
|
---|
945 | }
|
---|
946 |
|
---|
947 | /********* Methode *********/
|
---|
948 | //++
|
---|
949 | void Histo::HInteg(float norm)
|
---|
950 | //
|
---|
951 | // Remplace l'histogramme par son integrale normalise a norm:
|
---|
952 | // si norm <= 0 : pas de normalisation, integration seule
|
---|
953 | //--
|
---|
954 | {
|
---|
955 | if(bins>1)
|
---|
956 | for(int i=1; i<bins; i++) {
|
---|
957 | data[i] += data[i-1];
|
---|
958 | if(err2!=NULL) err2[i] += err2[i-1];
|
---|
959 | }
|
---|
960 | if( norm <= 0. ) return;
|
---|
961 | norm /= data[bins-1];
|
---|
962 | for(int i=0; i<bins; i++) {
|
---|
963 | data[i] *= norm;
|
---|
964 | if(err2!=NULL) err2[i] *= norm*norm;
|
---|
965 | }
|
---|
966 | }
|
---|
967 |
|
---|
968 | /********* Methode *********/
|
---|
969 | //++
|
---|
970 | void Histo::HDeriv()
|
---|
971 | //
|
---|
972 | // Remplace l'histogramme par sa derivee
|
---|
973 | //--
|
---|
974 | {
|
---|
975 | if( bins <= 1 ) return;
|
---|
976 | float *temp = new float[bins];
|
---|
977 | memcpy(temp, data, bins*sizeof(float));
|
---|
978 | if(bins>=3) for(int i=1; i<bins-1; i++)
|
---|
979 | temp[i] = (data[i+1]-data[i-1])/(2.*binWidth);
|
---|
980 | temp[0] = (data[1]-data[0])/binWidth;
|
---|
981 | temp[bins-1] = (data[bins-1]-data[bins-2])/binWidth;
|
---|
982 | memcpy(data, temp, bins*sizeof(float));
|
---|
983 | delete [] temp;
|
---|
984 | }
|
---|
985 |
|
---|
986 | /********* Methode *********/
|
---|
987 | //++
|
---|
988 | void Histo::HRebin(int nbinew)
|
---|
989 | //
|
---|
990 | // Pour rebinner l'histogramme sur nbinew bins
|
---|
991 | //--
|
---|
992 | {
|
---|
993 | if( nbinew <= 0 ) return;
|
---|
994 |
|
---|
995 | // memoraisation de l'histogramme original
|
---|
996 | Histo H(*this);
|
---|
997 |
|
---|
998 | // Le nombre de bins est il plus grand ??
|
---|
999 | if( nbinew > bins ) {
|
---|
1000 | delete [] data; data = NULL;
|
---|
1001 | data = new float[nbinew];
|
---|
1002 | if( err2 != NULL ) {
|
---|
1003 | delete [] err2; err2 = NULL;
|
---|
1004 | err2 = new double[nbinew];
|
---|
1005 | }
|
---|
1006 | }
|
---|
1007 |
|
---|
1008 | // remise en forme de this
|
---|
1009 | bins = nbinew;
|
---|
1010 | this->Zero();
|
---|
1011 | binWidth = (max-min)/bins;
|
---|
1012 | // On recopie les parametres qui n'ont pas changes
|
---|
1013 | under = H.under;
|
---|
1014 | over = H.over;
|
---|
1015 |
|
---|
1016 |
|
---|
1017 | // Remplissage
|
---|
1018 | for(int i=0;i<bins;i++) {
|
---|
1019 | float xmi = BinLowEdge(i);
|
---|
1020 | float xma = BinHighEdge(i);
|
---|
1021 | int imi = H.FindBin(xmi);
|
---|
1022 | if( imi < 0 ) imi = 0;
|
---|
1023 | int ima = H.FindBin(xma);
|
---|
1024 | if( ima >= H.bins ) ima = H.bins-1;
|
---|
1025 | double w = 0.;
|
---|
1026 | if( ima == imi ) {
|
---|
1027 | w = H(imi) * binWidth/H.binWidth;
|
---|
1028 | } else {
|
---|
1029 | w += H(imi) * (H.BinHighEdge(imi)-xmi)/H.binWidth;
|
---|
1030 | w += H(ima) * (xma -H.BinLowEdge(ima))/H.binWidth;
|
---|
1031 | if( ima > imi+1 )
|
---|
1032 | for(int ii=imi+1;ii<ima;ii++) w += H(ii);
|
---|
1033 | }
|
---|
1034 | AddBin(i,(float) w);
|
---|
1035 | }
|
---|
1036 |
|
---|
1037 | }
|
---|
1038 |
|
---|
1039 | /********* Methode *********/
|
---|
1040 | //++
|
---|
1041 | int Histo::MaxiLocal(float& maxi,int& imax,float& maxn,int& imaxn)
|
---|
1042 | //
|
---|
1043 | // Retourne le nombre de maxima locaux, la valeur du maximum (maxi)
|
---|
1044 | // et sa position (imax), ainsi que la valeur du second maximum
|
---|
1045 | // local (maxn) et sa position (imaxn).
|
---|
1046 | // Attention: si un maximum a la meme valeur sur plusieurs bins
|
---|
1047 | // consecutifs, le bin le plus a droite est pris.
|
---|
1048 | //--
|
---|
1049 | {
|
---|
1050 | int nml = 0;
|
---|
1051 | imax = imaxn = -1;
|
---|
1052 | maxi = maxn = -1.f;
|
---|
1053 |
|
---|
1054 | bool up = true;
|
---|
1055 | bool down = false;
|
---|
1056 | for(int i=0;i<bins;i++) {
|
---|
1057 | if( !up ) if( data[i] > data[i-1] ) up = true;
|
---|
1058 | if( up && !down ) {
|
---|
1059 | if( i == bins-1 ) down=true;
|
---|
1060 | else if( data[i] > data[i+1] ) down=true;
|
---|
1061 | }
|
---|
1062 |
|
---|
1063 | if( up && down ) {
|
---|
1064 | nml++;
|
---|
1065 | up = down = false;
|
---|
1066 | if( imax < 0 ) {
|
---|
1067 | imax = i; maxi = data[i];
|
---|
1068 | } else if( data[i] >= maxi ) {
|
---|
1069 | imaxn = imax; maxn = maxi;
|
---|
1070 | imax = i; maxi = data[i];
|
---|
1071 | } else {
|
---|
1072 | if( imaxn < 0 || data[i] >= maxn ) { imaxn = i; maxn = data[i]; }
|
---|
1073 | }
|
---|
1074 | }
|
---|
1075 |
|
---|
1076 | }
|
---|
1077 | return nml;
|
---|
1078 | }
|
---|
1079 |
|
---|
1080 | /********* Methode *********/
|
---|
1081 | //++
|
---|
1082 | float Histo::FitMax(int degree, float frac, int debug) const
|
---|
1083 | //
|
---|
1084 | // Fit de la position du maximum de l'histo par un polynome
|
---|
1085 | // de degre `degree' a `frac' pourcent du maximum.
|
---|
1086 | // L'histo est suppose etre remplit de valeurs positives
|
---|
1087 | //--
|
---|
1088 | {
|
---|
1089 | if (degree < 2 || degree > 3) degree = 3;
|
---|
1090 | if (frac <= 0. || frac >= 1.) frac = 0.5;
|
---|
1091 |
|
---|
1092 | if (debug > 1)
|
---|
1093 | cout<<"Histo::FitMax : Nb Entrees histo ="<<NEntries()<<endl;
|
---|
1094 |
|
---|
1095 | if (NEntries() < 1) THROW(inconsistentErr);
|
---|
1096 |
|
---|
1097 | int iMax = IMax();
|
---|
1098 | float hmax = (*this)(iMax);
|
---|
1099 | float xCenter = BinCenter(iMax);
|
---|
1100 |
|
---|
1101 | if(debug>1)
|
---|
1102 | cout<<"Max histo i="<<iMax<<" x="<<xCenter<<" v="<<hmax<<endl;
|
---|
1103 |
|
---|
1104 | /* find limits at frac*hmax */
|
---|
1105 |
|
---|
1106 | float limit = frac*hmax;
|
---|
1107 |
|
---|
1108 | volatile int iLow = iMax;
|
---|
1109 | while (iLow>0 && (*this)(iLow)>limit) iLow--;
|
---|
1110 |
|
---|
1111 | volatile int iHigh = iMax;
|
---|
1112 | while (iHigh<NBins()-1 && (*this)(iHigh)>limit) iHigh++;
|
---|
1113 |
|
---|
1114 | int nLowHigh;
|
---|
1115 | for(;;) {
|
---|
1116 | nLowHigh = 0;
|
---|
1117 | for (int i=iLow; i<=iHigh; i++) if((*this)(i)>0) {
|
---|
1118 | if(!err2) nLowHigh++;
|
---|
1119 | else if(Error2(i)>0.) nLowHigh++;
|
---|
1120 | }
|
---|
1121 | if (debug > 1) cout <<"Limites histo "<<iLow<<" - "<<iHigh
|
---|
1122 | <<" ("<<nLowHigh<<" non nuls)"<<endl;
|
---|
1123 | if( nLowHigh >= degree+1 ) break;
|
---|
1124 | iLow--; iHigh++;
|
---|
1125 | if(iLow<0 && iHigh>=NBins()) {
|
---|
1126 | if (debug>1)
|
---|
1127 | cout<<"Mode histogramme = "<<xCenter
|
---|
1128 | <<" BinCenter("<<iMax<<")"<<endl;
|
---|
1129 | return xCenter;
|
---|
1130 | }
|
---|
1131 | if(iLow < 0 ) iLow = 0;
|
---|
1132 | if(iHigh >= NBins()) iHigh = NBins()-1;
|
---|
1133 | }
|
---|
1134 |
|
---|
1135 | Vector xFit(nLowHigh);
|
---|
1136 | Vector yFit(nLowHigh);
|
---|
1137 | Vector e2Fit(nLowHigh);
|
---|
1138 | Vector errcoef(1);
|
---|
1139 | int ii = 0;
|
---|
1140 | for (int j=iLow; j<=iHigh; j++) {
|
---|
1141 | if ((*this)(j)>0) {
|
---|
1142 | if(!err2) {
|
---|
1143 | xFit(ii) = BinCenter(j)-xCenter;
|
---|
1144 | yFit(ii) = (*this)(j);
|
---|
1145 | e2Fit(ii) = yFit(ii);
|
---|
1146 | ii++;
|
---|
1147 | } else if(Error2(j)>0.) {
|
---|
1148 | xFit(ii) = BinCenter(j)-xCenter;
|
---|
1149 | yFit(ii) = (*this)(j);
|
---|
1150 | e2Fit(ii) = Error2(j);
|
---|
1151 | ii++;
|
---|
1152 | }
|
---|
1153 | }
|
---|
1154 | }
|
---|
1155 | if(debug>4) {
|
---|
1156 | int k;
|
---|
1157 | for(k=0;k<nLowHigh;k++) cout<<" "<<xFit(k); cout<<endl;
|
---|
1158 | for(k=0;k<nLowHigh;k++) cout<<" "<<yFit(k); cout<<endl;
|
---|
1159 | for(k=0;k<nLowHigh;k++) cout<<" "<<e2Fit(k); cout<<endl;
|
---|
1160 | }
|
---|
1161 | if( ii != nLowHigh ) THROW(inconsistentErr);
|
---|
1162 | Poly pol(degree);
|
---|
1163 | TRY {
|
---|
1164 | pol.Fit(xFit, yFit, e2Fit, degree, errcoef);
|
---|
1165 | if (debug>1) cout << "resultat fit = " << pol << endl;
|
---|
1166 | pol.Derivate();
|
---|
1167 | } CATCHALL {
|
---|
1168 | THROW_SAME;
|
---|
1169 | } ENDTRY
|
---|
1170 |
|
---|
1171 | int DPolDeg = pol.Degre();
|
---|
1172 | float fd = 0.;
|
---|
1173 |
|
---|
1174 | if (DPolDeg == 0) {
|
---|
1175 | // on est dans le cas d'un fit de droite
|
---|
1176 | if( pol[0] > 0. ) {
|
---|
1177 | // on a fitte une droite de pente >0.
|
---|
1178 | fd = xFit(nLowHigh-1) + binWidth/2. + xCenter;
|
---|
1179 | } else if( pol[0] < 0. ) {
|
---|
1180 | // on a fitte une droite de pente <0.
|
---|
1181 | fd = xFit(0) - binWidth/2. + xCenter;
|
---|
1182 | } else {
|
---|
1183 | // on a fitte une droite de pente =0. (creneau)
|
---|
1184 | fd = (xFit(0)+xFit(nLowHigh-1))/2. + xCenter;
|
---|
1185 | }
|
---|
1186 | } else if (DPolDeg == 1) {
|
---|
1187 | // on est dans le cas d'un fit de parabole
|
---|
1188 | double r=0;
|
---|
1189 | if (pol.Root1(r)==0) THROW(inconsistentErr);
|
---|
1190 | fd = r + xCenter;
|
---|
1191 | } else if (DPolDeg == 2) {
|
---|
1192 | // on est dans le cas d'un fit de cubique
|
---|
1193 | double r1=0;
|
---|
1194 | double r2=0;
|
---|
1195 | if (pol.Root2(r1,r2) == 0) THROW(inconsistentErr);
|
---|
1196 | pol.Derivate();
|
---|
1197 | fd = (pol(r1)<0) ? r1 + xCenter : r2 + xCenter;
|
---|
1198 | } else {
|
---|
1199 | // on est dans un cas non prevu
|
---|
1200 | THROW(inconsistentErr);
|
---|
1201 | }
|
---|
1202 |
|
---|
1203 | if(fd>max) fd = max;
|
---|
1204 | if(fd<min) fd = min;
|
---|
1205 |
|
---|
1206 | if (debug>1)
|
---|
1207 | cout << "Mode histogramme = " << fd
|
---|
1208 | << " (DerPol.degre =" << DPolDeg
|
---|
1209 | << " pol.coeff[0] =" << pol[0]
|
---|
1210 | << ")" << endl;
|
---|
1211 |
|
---|
1212 | return fd;
|
---|
1213 | }
|
---|
1214 |
|
---|
1215 |
|
---|
1216 | /********* Methode *********/
|
---|
1217 | //++
|
---|
1218 | float Histo::FindWidth(float frac, int debug) const
|
---|
1219 | //
|
---|
1220 | // Calcul de la largeur a frac pourcent du maximum
|
---|
1221 | // autour du bin du maximum.
|
---|
1222 | // L'histo est suppose etre remplit de valeurs positives
|
---|
1223 | //--
|
---|
1224 | {
|
---|
1225 | float xmax = BinCenter( IMax() );
|
---|
1226 | return FindWidth(xmax,frac,debug);
|
---|
1227 | }
|
---|
1228 |
|
---|
1229 | /********* Methode *********/
|
---|
1230 | //++
|
---|
1231 | float Histo::FindWidth(float xmax,float frac, int debug) const
|
---|
1232 | //
|
---|
1233 | // Calcul de la largeur a frac pourcent de la valeur du bin xmax.
|
---|
1234 | // L'histo est suppose etre remplit de valeurs positives
|
---|
1235 | //--
|
---|
1236 | {
|
---|
1237 | if (frac <= 0 || frac >= 1.) frac = 0.5;
|
---|
1238 |
|
---|
1239 | if (debug > 1)
|
---|
1240 | cout << "Histo::FindWidth a " << frac
|
---|
1241 | << " de xmax= " << xmax
|
---|
1242 | << " , ndata= " << NData()
|
---|
1243 | << " , nent= " << NEntries()
|
---|
1244 | << " , nbin= " << NBins() << endl;
|
---|
1245 |
|
---|
1246 | if (NEntries() < 1) THROW(inconsistentErr);
|
---|
1247 | if (NBins() < 3) THROW(inconsistentErr);
|
---|
1248 |
|
---|
1249 | int iMax = FindBin(xmax);
|
---|
1250 | if (iMax<0 || iMax>=NBins()) THROW(inconsistentErr);
|
---|
1251 | float hmax = data[iMax];
|
---|
1252 | float limit = frac*hmax;
|
---|
1253 | if (debug > 1)
|
---|
1254 | cout << " Max histo[" << iMax << "] = " << hmax
|
---|
1255 | << ", limite " << limit << endl;
|
---|
1256 |
|
---|
1257 | int iLow = iMax;
|
---|
1258 | while (iLow>=0 && data[iLow]>limit) iLow--;
|
---|
1259 | if( iLow < 0 ) iLow = 0;
|
---|
1260 |
|
---|
1261 | int iHigh = iMax;
|
---|
1262 | while (iHigh<NBins() && data[iHigh]>limit) iHigh++;
|
---|
1263 | if( iHigh >=NBins() ) iHigh = NBins()-1;
|
---|
1264 |
|
---|
1265 | float xlow = BinCenter(iLow);
|
---|
1266 | float ylow = data[iLow];
|
---|
1267 |
|
---|
1268 | float xlow1=xlow, ylow1=ylow;
|
---|
1269 | if(iLow+1<NBins()) {
|
---|
1270 | xlow1 = BinCenter(iLow+1);
|
---|
1271 | ylow1 = data[iLow+1];
|
---|
1272 | }
|
---|
1273 |
|
---|
1274 | float xhigh = BinCenter(iHigh);
|
---|
1275 | float yhigh = data[iHigh];
|
---|
1276 |
|
---|
1277 | float xhigh1=xhigh, yhigh1=yhigh;
|
---|
1278 | if(iHigh-1>=0) {
|
---|
1279 | xhigh1 = BinCenter(iHigh-1);
|
---|
1280 | yhigh1 = data[iHigh-1];
|
---|
1281 | }
|
---|
1282 |
|
---|
1283 | float xlpred,xhpred,wd;
|
---|
1284 |
|
---|
1285 | if(ylow1>ylow) {
|
---|
1286 | xlpred = xlow + (xlow1-xlow)/(ylow1-ylow)*(limit-ylow);
|
---|
1287 | if(xlpred < xlow) xlpred = xlow;
|
---|
1288 | } else xlpred = xlow;
|
---|
1289 |
|
---|
1290 | if(yhigh1>yhigh) {
|
---|
1291 | xhpred = xhigh + (xhigh1-xhigh)/(yhigh1-yhigh)*(limit-yhigh);
|
---|
1292 | if(xhpred > xhigh) xhpred = xhigh;
|
---|
1293 | } else xhpred = xhigh;
|
---|
1294 |
|
---|
1295 | wd = xhpred - xlpred;
|
---|
1296 |
|
---|
1297 | if (debug > 1) {
|
---|
1298 | cout << "Limites histo: " << " Width " << wd << endl;
|
---|
1299 | cout << " Low: [" << iLow << "]=" << ylow << " pred-> " << xlpred
|
---|
1300 | << " - High: [" << iHigh << "]=" << yhigh << " pred-> " << xhpred
|
---|
1301 | << endl;
|
---|
1302 | }
|
---|
1303 |
|
---|
1304 | return wd;
|
---|
1305 | }
|
---|
1306 |
|
---|
1307 |
|
---|
1308 | /********* Methode *********/
|
---|
1309 | //++
|
---|
1310 | int Histo::EstimeMax(float& xm,int SzPav)
|
---|
1311 | //
|
---|
1312 | // Cf suivant mais im est le bin du maximum de l'histo
|
---|
1313 | //--
|
---|
1314 | {
|
---|
1315 | int im = IMax();
|
---|
1316 | return EstimeMax(im,xm,SzPav);
|
---|
1317 | }
|
---|
1318 |
|
---|
1319 | /********* Methode *********/
|
---|
1320 | //++
|
---|
1321 | int Histo::EstimeMax(int& im,float& xm,int SzPav)
|
---|
1322 | //
|
---|
1323 | // Determine l'abscisses du maximum donne par im
|
---|
1324 | // en moyennant dans un pave SzPav autour du maximum
|
---|
1325 | //| Return:
|
---|
1326 | //| 0 = si fit maximum reussi avec SzPav pixels
|
---|
1327 | //| 1 = si fit maximum reussi avec moins que SzPav pixels
|
---|
1328 | //| 2 = si fit maximum echoue et renvoit BinCenter()
|
---|
1329 | //| -1 = si echec: SzPav <= 0 ou im hors limites
|
---|
1330 | //--
|
---|
1331 | {
|
---|
1332 | xm = 0;
|
---|
1333 | if( SzPav <= 0 ) return -1;
|
---|
1334 | if( im < 0 || im >= bins ) return -1;
|
---|
1335 |
|
---|
1336 | if( SzPav%2 == 0 ) SzPav++;
|
---|
1337 | SzPav = (SzPav-1)/2;
|
---|
1338 |
|
---|
1339 | int rc = 0;
|
---|
1340 | double dxm = 0, wx = 0;
|
---|
1341 | for(int i=im-SzPav;i<=im+SzPav;i++) {
|
---|
1342 | if( i<0 || i>= bins ) {rc=1; continue;}
|
---|
1343 | dxm += BinCenter(i) * (*this)(i);
|
---|
1344 | wx += (*this)(i);
|
---|
1345 | }
|
---|
1346 |
|
---|
1347 | if( wx > 0. ) {
|
---|
1348 | xm = dxm/wx;
|
---|
1349 | return rc;
|
---|
1350 | } else {
|
---|
1351 | xm = BinCenter(im);
|
---|
1352 | return 2;
|
---|
1353 | }
|
---|
1354 |
|
---|
1355 | }
|
---|
1356 |
|
---|
1357 | /********* Methode *********/
|
---|
1358 | //++
|
---|
1359 | void Histo::EstimeWidthS(float frac,float& widthG,float& widthD)
|
---|
1360 | //
|
---|
1361 | // Determine la largeur a frac% du maximum a gauche (widthG)
|
---|
1362 | // et a droite (widthD)
|
---|
1363 | //--
|
---|
1364 | {
|
---|
1365 | int i;
|
---|
1366 | widthG = widthD = -1.;
|
---|
1367 | if( bins<=1 || frac<=0. || frac>=1. ) return;
|
---|
1368 |
|
---|
1369 | int imax = 0;
|
---|
1370 | float maxi = data[0];
|
---|
1371 | for(i=1;i<bins;i++) if(data[i]>maxi) {imax=i; maxi=data[i];}
|
---|
1372 | float xmax = BinCenter(imax);
|
---|
1373 | float maxf = maxi * frac;
|
---|
1374 |
|
---|
1375 | // recherche du sigma a gauche
|
---|
1376 | widthG = 0.;
|
---|
1377 | for(i=imax;i>=0;i--) if( data[i] <= maxf ) break;
|
---|
1378 | if(i<0) i=0;
|
---|
1379 | if(i<imax) {
|
---|
1380 | if( data[i+1] != data[i] ) {
|
---|
1381 | widthG = BinCenter(i) + binWidth
|
---|
1382 | * (maxf-data[i])/(data[i+1]-data[i]);
|
---|
1383 | widthG = xmax - widthG;
|
---|
1384 | } else widthG = xmax - BinCenter(i);
|
---|
1385 | }
|
---|
1386 |
|
---|
1387 | // recherche du sigma a droite
|
---|
1388 | widthD = 0.;
|
---|
1389 | for(i=imax;i<bins;i++) if( data[i] <= maxf ) break;
|
---|
1390 | if(i>=bins) i=bins-1;
|
---|
1391 | if(i>imax) {
|
---|
1392 | if( data[i] != data[i-1] ) {
|
---|
1393 | widthD = BinCenter(i) - binWidth
|
---|
1394 | * (maxf-data[i])/(data[i-1]-data[i]);
|
---|
1395 | widthD -= xmax;
|
---|
1396 | } else widthD = BinCenter(i) - xmax;
|
---|
1397 | }
|
---|
1398 |
|
---|
1399 | }
|
---|
1400 |
|
---|
1401 | //////////////////////////////////////////////////////////
|
---|
1402 | //++
|
---|
1403 | int Histo::Fit(GeneralFit& gfit,unsigned short typ_err)
|
---|
1404 | //
|
---|
1405 | // Fit de l'histogramme par ``gfit''.
|
---|
1406 | //| typ_err = 0 :
|
---|
1407 | //| - erreur attachee au bin si elle existe
|
---|
1408 | //| - sinon 1
|
---|
1409 | //| typ_err = 1 :
|
---|
1410 | //| - erreur attachee au bin si elle existe
|
---|
1411 | //| - sinon max( sqrt(abs(bin) ,1 )
|
---|
1412 | //| typ_err = 2 :
|
---|
1413 | //| - erreur forcee a 1
|
---|
1414 | //| typ_err = 3 :
|
---|
1415 | //| - erreur forcee a max( sqrt(abs(bin) ,1 )
|
---|
1416 | //| typ_err = 4 :
|
---|
1417 | //| - erreur forcee a 1, nulle si bin a zero.
|
---|
1418 | //| typ_err = 5 :
|
---|
1419 | //| - erreur forcee a max( sqrt(abs(bin) ,1 ),
|
---|
1420 | //| nulle si bin a zero.
|
---|
1421 | //--
|
---|
1422 | {
|
---|
1423 | if(NBins()<=0) return -1000;
|
---|
1424 | if(typ_err>5) typ_err=0;
|
---|
1425 |
|
---|
1426 | GeneralFitData mydata(1,NBins());
|
---|
1427 |
|
---|
1428 | for(int i=0;i<NBins();i++) {
|
---|
1429 | double x = (double) BinCenter(i);
|
---|
1430 | double f = (double) (*this)(i);
|
---|
1431 | double saf = sqrt(fabs((double) f)); if(saf<1.) saf=1.;
|
---|
1432 | double e;
|
---|
1433 | if(typ_err==0) {if(HasErrors()) e=Error(i); else e=1.;}
|
---|
1434 | else if(typ_err==1) {if(HasErrors()) e=Error(i); else e=saf;}
|
---|
1435 | else if(typ_err==2) e=1.;
|
---|
1436 | else if(typ_err==3) e=saf;
|
---|
1437 | else if(typ_err==4) e=(f==0.)?0.:1.;
|
---|
1438 | else if(typ_err==5) e=(f==0.)?0.:saf;
|
---|
1439 | mydata.AddData1(x,f,e);
|
---|
1440 | }
|
---|
1441 |
|
---|
1442 | gfit.SetData(&mydata);
|
---|
1443 |
|
---|
1444 | return gfit.Fit();
|
---|
1445 | }
|
---|
1446 |
|
---|
1447 | //++
|
---|
1448 | Histo* Histo::FitResidus(GeneralFit& gfit)
|
---|
1449 | //
|
---|
1450 | // Retourne une classe contenant les residus du fit ``gfit''.
|
---|
1451 | //--
|
---|
1452 | {
|
---|
1453 | if(NBins()<=0) return NULL;
|
---|
1454 | GeneralFunction* f = gfit.GetFunction();
|
---|
1455 | if(f==NULL) return NULL;
|
---|
1456 | Vector par = gfit.GetParm();
|
---|
1457 | Histo* h = new Histo(*this);
|
---|
1458 | for(int i=0;i<NBins();i++) {
|
---|
1459 | double x = (double) BinCenter(i);
|
---|
1460 | (*h)(i) -= (float) f->Value(&x,par.Data());
|
---|
1461 | }
|
---|
1462 | return h;
|
---|
1463 | }
|
---|
1464 |
|
---|
1465 | //++
|
---|
1466 | Histo* Histo::FitFunction(GeneralFit& gfit)
|
---|
1467 | //
|
---|
1468 | // Retourne une classe contenant la fonction du fit ``gfit''.
|
---|
1469 | //--
|
---|
1470 | {
|
---|
1471 | if(NBins()<=0) return NULL;
|
---|
1472 | GeneralFunction* f = gfit.GetFunction();
|
---|
1473 | if(f==NULL) return NULL;
|
---|
1474 | Vector par = gfit.GetParm();
|
---|
1475 | Histo* h = new Histo(*this);
|
---|
1476 | for(int i=0;i<NBins();i++) {
|
---|
1477 | double x = (double) BinCenter(i);
|
---|
1478 | (*h)(i) = (float) f->Value(&x,par.Data());
|
---|
1479 | }
|
---|
1480 | return h;
|
---|
1481 | }
|
---|
1482 |
|
---|
1483 | /********* Methode *********/
|
---|
1484 | //++
|
---|
1485 | void Histo::PrintF(FILE * fp, int hdyn,float hmin, float hmax,int pflag,
|
---|
1486 | int il, int ih)
|
---|
1487 | //
|
---|
1488 | // Impression de l'histogramme dans le fichier fp
|
---|
1489 | //| hdyn = nombre de colonnes pour imprimer les etoiles
|
---|
1490 | //| si =0 alors defaut(100),
|
---|
1491 | //| si <0 pas de print histo seulement infos
|
---|
1492 | //| hmin = minimum de la dynamique
|
---|
1493 | //| hmax = maximum de la dynamique
|
---|
1494 | //| si hmax<=hmin : hmin=VMin() hmax=VMax()
|
---|
1495 | //| si hmax<=hmin et hmin=0 : hmin=0 hmax=VMax()
|
---|
1496 | //| sinon : hmin hmax
|
---|
1497 | //| pflag < 0 : on imprime les informations (nbin,min,...) sans l'histogramme
|
---|
1498 | //| = 0 : on imprime "BinCenter(i) data[i]" (note "... ...")
|
---|
1499 | //| bit 0 on : (v=1): numero_bin "... ..."
|
---|
1500 | //| bit 1 on : (v=2): "... ..." erreur
|
---|
1501 | //--
|
---|
1502 | {
|
---|
1503 |
|
---|
1504 |
|
---|
1505 | if( il > ih ) { il = 0; ih = bins-1; }
|
---|
1506 | if( il < 0 ) il = 0;
|
---|
1507 | if( ih >= bins ) ih = bins-1;
|
---|
1508 |
|
---|
1509 | double dhmin = (double) hmin;
|
---|
1510 | double dhmax = (double) hmax;
|
---|
1511 | double hb,hbmn,hbmx;
|
---|
1512 |
|
---|
1513 | if(hdyn==0) hdyn = 100;
|
---|
1514 |
|
---|
1515 | cout << "~Histo::Print "
|
---|
1516 | << " nHist=" << nHist << " nEntries=" << nEntries
|
---|
1517 | << " under=" << under << " over=" << over << endl;
|
---|
1518 | cout << " bins=" << bins
|
---|
1519 | << " min=" << min << " max=" << max
|
---|
1520 | << " binWidth=" << binWidth << endl;
|
---|
1521 | cout << " mean=" << Mean() << " r.m.s=" << Sigma() << endl;
|
---|
1522 |
|
---|
1523 | if(hdyn<0 || pflag<0 ) return;
|
---|
1524 |
|
---|
1525 | if(dhmax<=dhmin) { if(hmin != 0.) dhmin = (double) VMin(); else dhmin=0.;
|
---|
1526 | dhmax = (double) VMax(); }
|
---|
1527 | if(dhmin>dhmax) return;
|
---|
1528 | if(dhmin==dhmax) {dhmin -= 1.; dhmax += 1.;}
|
---|
1529 | double wb = (dhmax-dhmin) / (double) hdyn;
|
---|
1530 |
|
---|
1531 | // determination de la position de la valeur zero
|
---|
1532 | int i0 = (int)(-dhmin/wb);
|
---|
1533 |
|
---|
1534 | // si le zero est dans la dynamique,
|
---|
1535 | // il doit imperativement etre une limite de bin
|
---|
1536 | if( 0 <= i0 && i0 < hdyn ) {
|
---|
1537 | hbmn = dhmin + i0*wb;
|
---|
1538 | hbmx = hbmn + wb;
|
---|
1539 | if( hbmn != 0. ) {
|
---|
1540 | hbmn *= -1.;
|
---|
1541 | if( hbmn < hbmx ) {
|
---|
1542 | // le zero est + pres du bord negatif du bin
|
---|
1543 | dhmin += hbmn;
|
---|
1544 | dhmax += hbmn;
|
---|
1545 | } else {
|
---|
1546 | // le zero est + pres du bord positif du bin
|
---|
1547 | dhmin -= hbmx;
|
---|
1548 | dhmax -= hbmx;
|
---|
1549 | }
|
---|
1550 | wb = (dhmax-dhmin) / (double) hdyn;
|
---|
1551 | i0 = (int)(-dhmin/wb);
|
---|
1552 | }
|
---|
1553 | }
|
---|
1554 |
|
---|
1555 | cout <<" bin minimum="<<dhmin<<" bin maximum="<<dhmax
|
---|
1556 | <<" binw="<<wb<< endl;
|
---|
1557 |
|
---|
1558 | char* s = new char[hdyn+1];
|
---|
1559 | s[hdyn] = '\0';
|
---|
1560 |
|
---|
1561 | // premiere ligne
|
---|
1562 | {for(int i=0;i<hdyn;i++) s[i] = '=';}
|
---|
1563 | if( 0 <= i0 && i0 < hdyn ) s[i0] = '0';
|
---|
1564 | if(pflag&1) fprintf( fp, "====");
|
---|
1565 | fprintf( fp, "======================");
|
---|
1566 | if(pflag&2 && err2!=NULL) fprintf( fp, "===========");
|
---|
1567 | fprintf( fp, "==%s\n",s);
|
---|
1568 |
|
---|
1569 | // histogramme
|
---|
1570 | {for(int i=il;i<=ih;i++) {
|
---|
1571 | for(int k=0;k<hdyn;k++) s[k] = ' ';
|
---|
1572 | hb = (*this)(i);
|
---|
1573 |
|
---|
1574 | //choix du bin (hdyn bin entre [dhmin,dhmax[
|
---|
1575 | int ii = (int)( (hb-dhmin)/wb );
|
---|
1576 | if(ii<0) ii = 0; else if(ii>=hdyn) ii = hdyn-1;
|
---|
1577 |
|
---|
1578 | // limite du bin
|
---|
1579 | hbmn = dhmin + ii*wb;
|
---|
1580 | hbmx = hbmn + wb;
|
---|
1581 |
|
---|
1582 | // remplissage de s[] en tenant compte des courbes +/-
|
---|
1583 | if(i0<0) {
|
---|
1584 | // courbe entierement positive
|
---|
1585 | for(int k=0;k<=ii;k++) s[k] = 'X';
|
---|
1586 | } else if(i0>=hdyn) {
|
---|
1587 | // courbe entierement negative
|
---|
1588 | for(int k=hdyn-1;k>=ii;k--) s[k] = 'X';
|
---|
1589 | } else {
|
---|
1590 | // courbe positive et negative
|
---|
1591 | s[i0] = '|';
|
---|
1592 | if(ii>i0) for(int k=i0+1;k<=ii;k++) s[k] = 'X';
|
---|
1593 | else if(ii<i0) for(int k=ii;k<i0;k++) s[k] = 'X';
|
---|
1594 | else s[ii] = 'X';
|
---|
1595 | }
|
---|
1596 |
|
---|
1597 | // valeur a mettre dans le bin le plus haut/bas
|
---|
1598 | int ib;
|
---|
1599 | if(hb>0.) ib = (int)( 10.*(hb-hbmn)/(hbmx-hbmn) );
|
---|
1600 | else if(hb<0.) ib = (int)( 10.*(hbmx-hb)/(hbmx-hbmn) );
|
---|
1601 | else ib = -1;
|
---|
1602 | if(ib==-1) s[ii] = '|';
|
---|
1603 | else if(ib==0) s[ii] = '.';
|
---|
1604 | else if(ib>0 && ib<10) s[ii] = (char)((int) '0' + ib);
|
---|
1605 |
|
---|
1606 | // traitement des saturations +/-
|
---|
1607 | if( hb < dhmin ) s[0] = '*';
|
---|
1608 | else if( hb > dhmax ) s[hdyn-1] = '*';
|
---|
1609 |
|
---|
1610 | if(pflag&1) fprintf( fp, "%3d ",i);
|
---|
1611 | fprintf( fp, "%10.4g %10.4g ",BinCenter(i),hb);
|
---|
1612 | if(pflag&2 && err2!=NULL) fprintf( fp, "%10.4g ",Error(i));
|
---|
1613 | fprintf( fp, "= %s\n",s);
|
---|
1614 | }}
|
---|
1615 |
|
---|
1616 | // derniere ligne
|
---|
1617 | for(int i=0;i<hdyn;i++) s[i] = '=';
|
---|
1618 | if( 0 <= i0 && i0 < hdyn ) s[i0] = '0';
|
---|
1619 | if(pflag&1) fprintf( fp, "====");
|
---|
1620 | fprintf( fp, "======================");
|
---|
1621 | if(pflag&2 && err2!=NULL) fprintf( fp, "===========");
|
---|
1622 | fprintf( fp, "==%s\n",s);
|
---|
1623 |
|
---|
1624 | // valeur basse des bins (sur ["ndig-1" digits + signe] = ndig char (>=3))
|
---|
1625 | const int ndig = 7;
|
---|
1626 | char sn[2*ndig+10];
|
---|
1627 | hb = ( fabs(dhmin) > fabs(dhmax) ) ? fabs(dhmin) : fabs(dhmax);
|
---|
1628 | int n;
|
---|
1629 | if( hb > 0. ) n = (int)(log10(hb*1.00000001)); else n = 1;
|
---|
1630 | double scaledig = pow(10.,(double) ndig-2);
|
---|
1631 | double expo = scaledig/pow(10.,(double) n);
|
---|
1632 | // cout <<"n="<<n<<" hb="<<hb<<" scaledig="<<scaledig<<" expo="<<expo<<endl;
|
---|
1633 | for(int k=0;k<ndig;k++) {
|
---|
1634 | for(int i=0;i<hdyn;i++) s[i] = ' ';
|
---|
1635 | {for(int i=0;i<hdyn;i++) {
|
---|
1636 | n = (int)( (dhmin + i*wb)*expo );
|
---|
1637 | for(int j=0;j<2*ndig+10;j++) sn[j] = ' ';
|
---|
1638 | sprintf(sn,"%d%c",n,'\0');
|
---|
1639 | strip(sn,'B',' ');
|
---|
1640 | // cout <<"n="<<n<<" sn=("<<sn<<") l="<<strlen(sn)<<" k="<<k;
|
---|
1641 | if( (int) strlen(sn) > k ) s[i] = sn[k];
|
---|
1642 | // cout <<" s=("<<s<<")"<<endl;
|
---|
1643 | }}
|
---|
1644 | if(pflag&1) fprintf( fp, " ");
|
---|
1645 | fprintf( fp, " ");
|
---|
1646 | if(pflag&2 && err2!=NULL) fprintf( fp, " ");
|
---|
1647 | fprintf( fp, " %s\n",s);
|
---|
1648 | }
|
---|
1649 | fprintf( fp, " (valeurs a multiplier par %.0e)\n",1./expo);
|
---|
1650 |
|
---|
1651 | delete[] s;
|
---|
1652 | }
|
---|
1653 |
|
---|
1654 | /********* Methode *********/
|
---|
1655 | //++
|
---|
1656 | void Histo::Print(int hdyn,float hmin, float hmax,int pflag,
|
---|
1657 | int il, int ih)
|
---|
1658 | //
|
---|
1659 | // Impression de l'histogramme sur stdout
|
---|
1660 | //--
|
---|
1661 | {
|
---|
1662 | Histo::PrintF(stdout, hdyn, hmin, hmax, pflag, il, ih);
|
---|
1663 | }
|
---|
1664 |
|
---|
1665 |
|
---|
1666 | /********* Methode *********/
|
---|
1667 | //++
|
---|
1668 | void Histo::WriteSelf(POutPersist& s) const
|
---|
1669 | //
|
---|
1670 | // Ecriture ppersist
|
---|
1671 | //--
|
---|
1672 | {
|
---|
1673 | char strg[256];
|
---|
1674 |
|
---|
1675 | // Que faut-il ecrire?
|
---|
1676 | int errok = (err2) ? 1 : 0;
|
---|
1677 |
|
---|
1678 | // Ecriture entete pour identifier facilement
|
---|
1679 | sprintf(strg,"bins=%6d NEnt=%15d errok=%1d",bins,nEntries,errok);
|
---|
1680 | s.PutLine(strg);
|
---|
1681 | sprintf(strg,"binw=%g min=%g max=%g",binWidth,min,max);
|
---|
1682 | s.PutLine(strg);
|
---|
1683 | sprintf(strg, "under=%g over=%g nHist=%g",under,over,nHist);
|
---|
1684 | s.PutLine(strg);
|
---|
1685 |
|
---|
1686 | // Ecriture des valeurs
|
---|
1687 | s.PutI4(bins);
|
---|
1688 | s.PutI4(nEntries);
|
---|
1689 | s.PutI4(errok);
|
---|
1690 |
|
---|
1691 | s.PutR4(binWidth);
|
---|
1692 | s.PutR4(min);
|
---|
1693 | s.PutR4(max);
|
---|
1694 | s.PutR4(under);
|
---|
1695 | s.PutR4(over);
|
---|
1696 |
|
---|
1697 | s.PutR8(nHist);
|
---|
1698 |
|
---|
1699 | // Ecriture des donnees Histo 1D
|
---|
1700 | sprintf(strg,"Histo: Tableau des donnees %d",bins);
|
---|
1701 | s.PutLine(strg);
|
---|
1702 | s.PutR4s(data, bins);
|
---|
1703 |
|
---|
1704 | // Ecriture des erreurs
|
---|
1705 | if(errok) {
|
---|
1706 | sprintf(strg,"Histo: Tableau des erreurs %d",bins);
|
---|
1707 | s.PutLine(strg);
|
---|
1708 | s.PutR8s(err2, bins);
|
---|
1709 | }
|
---|
1710 |
|
---|
1711 | return;
|
---|
1712 | }
|
---|
1713 |
|
---|
1714 | /********* Methode *********/
|
---|
1715 | //++
|
---|
1716 | void Histo::ReadSelf(PInPersist& s)
|
---|
1717 | //
|
---|
1718 | // Lecture ppersist
|
---|
1719 | //--
|
---|
1720 | {
|
---|
1721 | char strg[256];
|
---|
1722 |
|
---|
1723 | Delete();
|
---|
1724 |
|
---|
1725 | // Lecture entete
|
---|
1726 | s.GetLine(strg, 255);
|
---|
1727 | s.GetLine(strg, 255);
|
---|
1728 | s.GetLine(strg, 255);
|
---|
1729 |
|
---|
1730 | // Lecture des valeurs
|
---|
1731 | s.GetI4(bins);
|
---|
1732 | s.GetI4(nEntries);
|
---|
1733 | int_4 errok;
|
---|
1734 | s.GetI4(errok);
|
---|
1735 |
|
---|
1736 | s.GetR4(binWidth);
|
---|
1737 | s.GetR4(min);
|
---|
1738 | s.GetR4(max);
|
---|
1739 | s.GetR4(under);
|
---|
1740 | s.GetR4(over);
|
---|
1741 |
|
---|
1742 | s.GetR8(nHist);
|
---|
1743 |
|
---|
1744 | // Lecture des donnees Histo 1D
|
---|
1745 | data = new float[bins];
|
---|
1746 | s.GetLine(strg, 255);
|
---|
1747 | s.GetR4s(data, bins);
|
---|
1748 |
|
---|
1749 | // Lecture des erreurs
|
---|
1750 | if(errok) {
|
---|
1751 | s.GetLine(strg, 255);
|
---|
1752 | err2 = new double[bins];
|
---|
1753 | s.GetR8s(err2, bins);
|
---|
1754 | }
|
---|
1755 |
|
---|
1756 | return;
|
---|
1757 | }
|
---|
1758 |
|
---|
1759 |
|
---|
1760 |
|
---|
1761 | // Rappel des inlines functions pour commentaires
|
---|
1762 | //++
|
---|
1763 | // inline float XMin() const
|
---|
1764 | // Retourne l'abscisse minimum
|
---|
1765 | //--
|
---|
1766 | //++
|
---|
1767 | // inline float XMax() const
|
---|
1768 | // Retourne l'abscisse maximum
|
---|
1769 | //--
|
---|
1770 | //++
|
---|
1771 | // inline int NBins() const
|
---|
1772 | // Retourne le nombre de bins
|
---|
1773 | //--
|
---|
1774 | //++
|
---|
1775 | // inline float BinWidth() const
|
---|
1776 | // Retourne la largeur du bin
|
---|
1777 | //--
|
---|
1778 | //++
|
---|
1779 | // inline float* Bins() const
|
---|
1780 | // Retourne le pointeur sur le tableaux des contenus
|
---|
1781 | //--
|
---|
1782 | //++
|
---|
1783 | // inline float operator()(int i) const
|
---|
1784 | // Retourne le contenu du bin i
|
---|
1785 | //--
|
---|
1786 | //++
|
---|
1787 | // inline float& operator()(int i)
|
---|
1788 | // Remplit le contenu du bin i
|
---|
1789 | //--
|
---|
1790 | //++
|
---|
1791 | // inline float Error(int i) const
|
---|
1792 | // Retourne l'erreur du bin i
|
---|
1793 | //--
|
---|
1794 | //++
|
---|
1795 | // inline double Error2(int i) const
|
---|
1796 | // Retourne l'erreur au carre du bin i
|
---|
1797 | //--
|
---|
1798 | //++
|
---|
1799 | // inline double& Error2(int i)
|
---|
1800 | // Remplit l'erreur au carre du bin i
|
---|
1801 | //--
|
---|
1802 | //++
|
---|
1803 | // inline float NData() const
|
---|
1804 | // Retourne la somme ponderee
|
---|
1805 | //--
|
---|
1806 | //++
|
---|
1807 | // inline float NEntries()
|
---|
1808 | // Retourne le nombre d'entrees
|
---|
1809 | //--
|
---|
1810 | //++
|
---|
1811 | // inline float NOver() const
|
---|
1812 | // Retourne le nombre d'overflow
|
---|
1813 | //--
|
---|
1814 | //++
|
---|
1815 | // inline float NUnder() const
|
---|
1816 | // Retourne le nombre d'underflow
|
---|
1817 | //--
|
---|
1818 | //++
|
---|
1819 | // inline float BinLowEdge(int i) const
|
---|
1820 | // Retourne l'abscisse du bord inferieur du bin i
|
---|
1821 | //--
|
---|
1822 | //++
|
---|
1823 | // inline float BinCenter(int i) const
|
---|
1824 | // Retourne l'abscisse du centre du bin i
|
---|
1825 | //--
|
---|
1826 | //++
|
---|
1827 | // inline float BinHighEdge(int i) const
|
---|
1828 | // Retourne l'abscisse du bord superieur du bin i
|
---|
1829 | //--
|
---|
1830 | //++
|
---|
1831 | // inline int FindBin(float x) const
|
---|
1832 | // Retourne le numero du bin contenant l'abscisse x
|
---|
1833 | //--
|
---|