source: Sophya/trunk/SophyaLib/NTools/histos.cc@ 220

Last change on this file since 220 was 220, checked in by ansari, 26 years ago

Creation module DPC/NTools Reza 09/04/99

File size: 38.2 KB
Line 
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//++
30Histo::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//++
43Histo::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//++
57Histo::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//++
71Histo::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 *********/
88void 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 *********/
98Histo::~Histo()
99{
100Delete();
101}
102
103//++
104// Titre Methodes
105//--
106
107/********* Methode *********/
108//++
109void 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//++
123void 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//++
136Histo& 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//++
163Histo& Histo::operator *= (double b)
164//
165//--
166{
167double b2 = b*b;
168for(int i=0;i<bins;i++) {
169 data[i] *= b;
170 if(err2) err2[i] *= b2;
171}
172under *= b;
173over *= b;
174nHist *= b;
175
176return *this;
177}
178
179//++
180Histo& Histo::operator /= (double b)
181//
182//--
183{
184if (b==0.) THROW(inconsistentErr);
185double b2 = b*b;
186for(int i=0;i<bins;i++) {
187 data[i] /= b;
188 if(err2) err2[i] /= b2;
189}
190under /= b;
191over /= b;
192nHist /= b;
193
194return *this;
195}
196
197//++
198Histo& Histo::operator += (double b)
199//
200//--
201{
202for(int i=0;i<bins;i++) data[i] += b;
203under += b;
204over += b;
205nHist += bins*b;
206
207return *this;
208}
209
210//++
211Histo& Histo::operator -= (double b)
212//
213//--
214{
215for(int i=0;i<bins;i++) data[i] -= b;
216under -= b;
217over -= b;
218nHist -= bins*b;
219
220return *this;
221}
222
223/********* Methode *********/
224//++
225Histo operator * (const Histo& a, double b)
226//
227//--
228{
229 Histo result(a);
230 return (result *= b);
231}
232
233//++
234Histo operator * (double b, const Histo& a)
235//
236//--
237{
238 Histo result(a);
239 return (result *= b);
240}
241
242//++
243Histo operator / (const Histo& a, double b)
244//
245//--
246{
247 Histo result(a);
248 return (result /= b);
249}
250
251//++
252Histo operator + (const Histo& a, double b)
253//
254//--
255{
256 Histo result(a);
257 return (result += b);
258}
259
260//++
261Histo operator + (double b, const Histo& a)
262//
263//--
264{
265 Histo result(a);
266 return (result += b);
267}
268
269//++
270Histo operator - (const Histo& a, double b)
271//
272//--
273{
274 Histo result(a);
275 return (result -= b);
276}
277
278//++
279Histo 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//++
290Histo& Histo::operator += (const Histo& a)
291//
292//--
293{
294if(bins!=a.bins) THROW(sizeMismatchErr);
295for(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}
300under += a.under;
301over += a.over;
302nHist += a.nHist;
303nEntries += a.nEntries;
304
305return *this;
306}
307
308//++
309Histo& Histo::operator -= (const Histo& a)
310//
311//--
312{
313if(bins!=a.bins) THROW(sizeMismatchErr);
314for(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}
319under -= a.under;
320over -= a.over;
321nHist -= a.nHist;
322nEntries += a.nEntries;
323
324return *this;
325}
326
327//++
328Histo& Histo::operator *= (const Histo& a)
329//
330//--
331{
332if(bins!=a.bins) THROW(sizeMismatchErr);
333nHist = 0.;
334for(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}
341under *= a.under;
342over *= a.over;
343nEntries += a.nEntries;
344
345return *this;
346}
347
348//++
349Histo& Histo::operator /= (const Histo& a)
350//
351//--
352{
353if(bins!=a.bins) THROW(sizeMismatchErr);
354nHist = 0.;
355for(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}
368if(a.under!=0.) under /= a.under; else under = 0.;
369if(a.over!=0.) over /= a.over; else over = 0.;
370nEntries += a.nEntries;
371
372return *this;
373}
374
375/********* Methode *********/
376//++
377Histo operator + (const Histo& a, const Histo& b)
378//
379//--
380{
381if (b.bins!=a.bins) THROW(sizeMismatchErr);
382Histo c(a);
383return (c += b);
384}
385
386//++
387Histo operator - (const Histo& a, const Histo& b)
388//
389//--
390{
391if (b.bins!=a.bins) THROW(sizeMismatchErr);
392Histo c(a);
393return (c -= b);
394}
395
396//++
397Histo operator * (const Histo& a, const Histo& b)
398//
399//--
400{
401if (b.bins!=a.bins) THROW(sizeMismatchErr);
402Histo c(a);
403return (c *= b);
404}
405
406//++
407Histo operator / (const Histo& a, const Histo& b)
408//
409//--
410{
411if (b.bins!=a.bins) THROW(sizeMismatchErr);
412Histo c(a);
413return (c /= b);
414}
415
416/********* Methode *********/
417//++
418void Histo::GetAbsc(Vector &v)
419//
420// Remplissage d'un tableau avec la valeur des abscisses
421//--
422{
423v.Realloc(bins);
424for(int i=0;i<bins;i++) v(i) = BinLowEdge(i);
425return;
426}
427
428//++
429void Histo::GetValue(Vector &v)
430//
431// Remplissage d'un tableau avec la valeur du contenu
432//--
433{
434v.Realloc(bins);
435for(int i=0;i<bins;i++) v(i) = data[i];
436return;
437}
438
439//++
440void Histo::GetError2(Vector &v)
441//
442// Remplissage d'un tableau avec la valeur des erreurs au carre
443//--
444{
445v.Realloc(bins);
446if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
447for(int i=0;i<bins;i++) v(i) = err2[i];
448return;
449}
450
451//++
452void Histo::GetError(Vector &v)
453//
454// Remplissage d'un tableau avec la valeur des erreurs
455//--
456{
457v.Realloc(bins);
458if(!err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}
459for(int i=0;i<bins;i++) v(i) = Error(i);
460return;
461}
462
463/********* Methode *********/
464//++
465void Histo::PutValue(Vector &v, int ierr)
466//
467// Remplissage du contenu de l'histo avec les valeurs d'un vecteur
468//--
469{
470if(v.NElts()<bins) THROW(sizeMismatchErr);
471for(int i=0;i<bins;i++) {
472 data[i] = v(i);
473 if(err2&&ierr) err2[i] = fabs(v(i));
474}
475return;
476}
477
478//++
479void Histo::PutValueAdd(Vector &v, int ierr)
480//
481// Addition du contenu de l'histo avec les valeurs d'un vecteur
482//--
483{
484if(v.NElts()<bins) THROW(sizeMismatchErr);
485for(int i=0;i<bins;i++) {
486 data[i] += v(i);
487 if(err2&&ierr) err2[i] += fabs(v(i));
488}
489return;
490}
491
492//++
493void Histo::PutError2(Vector &v)
494//
495// Remplissage des erreurs au carre de l'histo avec les valeurs d'un vecteur
496//--
497{
498if(v.NElts()<bins) THROW(sizeMismatchErr);
499if(!err2) Errors();
500for(int i=0;i<bins;i++) err2[i] = v(i);
501return;
502}
503
504//++
505void Histo::PutError2Add(Vector &v)
506//
507// Addition des erreurs au carre de l'histo avec les valeurs d'un vecteur
508//--
509{
510if(v.NElts()<bins) THROW(sizeMismatchErr);
511if(!err2) Errors();
512for(int i=0;i<bins;i++) if(v(i)>0.) err2[i] += v(i);
513return;
514}
515
516//++
517void Histo::PutError(Vector &v)
518//
519// Remplissage des erreurs de l'histo avec les valeurs d'un vecteur
520//--
521{
522if(v.NElts()<bins) THROW(sizeMismatchErr);
523if(!err2) Errors();
524for(int i=0;i<bins;i++)
525 if(v(i)>0.) err2[i]=v(i)*v(i); else err2[i]=-v(i)*v(i);
526return;
527}
528
529/********* Methode *********/
530//++
531void 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//++
549void 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//++
566void 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//++
577void 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//++
593void 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//++
604void 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//++
616void Histo::SetErr(float x, float e)
617//
618// Remplissage des erreurs pour abscisse x
619//--
620{
621SetErr2(x, (double) e*e);
622}
623
624/********* Methode *********/
625//++
626void Histo::SetErr(int numBin, float e)
627//
628// Remplissage des erreurs pour numBin
629//--
630{
631SetErr2(numBin, (double) e*e);
632}
633
634/********* Methode *********/
635//++
636int 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//++
651int 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//++
666float 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//++
679float 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//++
692float 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//++
709float 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//++
732float 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//++
752float 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//++
777float 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//++
798float 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//++
822float 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//++
843float 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//++
855int Histo::BinNonNul() const
856//
857// Retourne le nombre de bins non-nul
858//--
859{
860int non=0;
861for (int i=0;i<bins;i++) if( data[i] != 0. ) non++;
862return non;
863}
864
865/********* Methode *********/
866//++
867int Histo::ErrNonNul() const
868//
869// Retourne le nombre de bins ayant une erreur non-nulle
870//--
871{
872if(err2==NULL) return -1;
873int non=0;
874for (int i=0;i<bins;i++) if( err2[i] != 0. ) non++;
875return non;
876}
877
878/********* Methode *********/
879//++
880int 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{
886double n = nHist*per;
887double s = 0.;
888int i;
889for(i=0; i<bins; i++ ) {
890 s += data[i];
891 if( s >= n ) break;
892}
893if( i == bins ) i = bins-1;
894return i;
895}
896
897/********* Methode *********/
898//++
899int 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{
911imin = imax = -1;
912if( per <= 0. ) return -1;
913
914int I = FindBin(x);
915if( I<0 || I>=bins ) return -2;
916
917double N1 = 0.; for(int i=0; i<=I; i++) N1 += data[i]; N1 *= per;
918double N2 = 0.; {for(int i=I; i<bins; i++) N2 += data[i]; N2 *= per;}
919
920double n = 0.;
921for(imin=I; imin>=0; imin--) { n += data[imin]; if(n>=N1) break; }
922if( imin<0 ) imin = 0;
923// cout<<"I="<<I<<" imin="<<imin<<" n="<<n<<" N1="<<N1<<endl;
924
925n = 0.;
926for(imax=I; imax<bins; imax++) { n += data[imax]; if(n>=N2) break; }
927if( imax>=bins ) imax = bins-1;
928// cout<<"I="<<I<<" imax="<<imax<<" n="<<n<<" N2="<<N2<<endl;
929
930return ( imax-I > I-imin ) ? I-imin : imax-I ;
931}
932
933/********* Methode *********/
934//++
935int Histo::BinPercent(float x,float per,float& xmin,float& xmax)
936//
937// Idem precedent mais renvoie xmin et xmax
938//--
939{
940xmin = xmax = 0.;
941int imin,imax;
942int rc = BinPercent(x,per,imin,imax);
943if( rc >= 0 ) { xmin = BinLowEdge(imin); xmax = BinHighEdge(imax); }
944return rc;
945}
946
947/********* Methode *********/
948//++
949void 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{
955if(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 }
960if( norm <= 0. ) return;
961norm /= data[bins-1];
962for(int i=0; i<bins; i++) {
963 data[i] *= norm;
964 if(err2!=NULL) err2[i] *= norm*norm;
965}
966}
967
968/********* Methode *********/
969//++
970void Histo::HDeriv()
971//
972// Remplace l'histogramme par sa derivee
973//--
974{
975if( bins <= 1 ) return;
976float *temp = new float[bins];
977memcpy(temp, data, bins*sizeof(float));
978if(bins>=3) for(int i=1; i<bins-1; i++)
979 temp[i] = (data[i+1]-data[i-1])/(2.*binWidth);
980temp[0] = (data[1]-data[0])/binWidth;
981temp[bins-1] = (data[bins-1]-data[bins-2])/binWidth;
982memcpy(data, temp, bins*sizeof(float));
983delete [] temp;
984}
985
986/********* Methode *********/
987//++
988void 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//++
1041int 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{
1050int nml = 0;
1051imax = imaxn = -1;
1052maxi = maxn = -1.f;
1053
1054bool up = true;
1055bool down = false;
1056for(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}
1077return nml;
1078}
1079
1080/********* Methode *********/
1081//++
1082float 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//++
1218float 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//++
1231float 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//++
1310int Histo::EstimeMax(float& xm,int SzPav)
1311//
1312// Cf suivant mais im est le bin du maximum de l'histo
1313//--
1314{
1315int im = IMax();
1316return EstimeMax(im,xm,SzPav);
1317}
1318
1319/********* Methode *********/
1320//++
1321int 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{
1332xm = 0;
1333if( SzPav <= 0 ) return -1;
1334if( im < 0 || im >= bins ) return -1;
1335
1336if( SzPav%2 == 0 ) SzPav++;
1337SzPav = (SzPav-1)/2;
1338
1339int rc = 0;
1340double dxm = 0, wx = 0;
1341for(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
1347if( 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//++
1359void 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{
1365int i;
1366widthG = widthD = -1.;
1367if( bins<=1 || frac<=0. || frac>=1. ) return;
1368
1369int imax = 0;
1370float maxi = data[0];
1371for(i=1;i<bins;i++) if(data[i]>maxi) {imax=i; maxi=data[i];}
1372float xmax = BinCenter(imax);
1373float maxf = maxi * frac;
1374
1375// recherche du sigma a gauche
1376widthG = 0.;
1377for(i=imax;i>=0;i--) if( data[i] <= maxf ) break;
1378if(i<0) i=0;
1379if(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
1388widthD = 0.;
1389for(i=imax;i<bins;i++) if( data[i] <= maxf ) break;
1390if(i>=bins) i=bins-1;
1391if(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//++
1403int 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{
1423if(NBins()<=0) return -1000;
1424if(typ_err>5) typ_err=0;
1425
1426GeneralFitData mydata(1,NBins());
1427
1428for(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
1442gfit.SetData(&mydata);
1443
1444return gfit.Fit();
1445}
1446
1447//++
1448Histo* Histo::FitResidus(GeneralFit& gfit)
1449//
1450// Retourne une classe contenant les residus du fit ``gfit''.
1451//--
1452{
1453if(NBins()<=0) return NULL;
1454GeneralFunction* f = gfit.GetFunction();
1455if(f==NULL) return NULL;
1456Vector par = gfit.GetParm();
1457Histo* h = new Histo(*this);
1458for(int i=0;i<NBins();i++) {
1459 double x = (double) BinCenter(i);
1460 (*h)(i) -= (float) f->Value(&x,par.Data());
1461}
1462return h;
1463}
1464
1465//++
1466Histo* Histo::FitFunction(GeneralFit& gfit)
1467//
1468// Retourne une classe contenant la fonction du fit ``gfit''.
1469//--
1470{
1471if(NBins()<=0) return NULL;
1472GeneralFunction* f = gfit.GetFunction();
1473if(f==NULL) return NULL;
1474Vector par = gfit.GetParm();
1475Histo* h = new Histo(*this);
1476for(int i=0;i<NBins();i++) {
1477 double x = (double) BinCenter(i);
1478 (*h)(i) = (float) f->Value(&x,par.Data());
1479}
1480return h;
1481}
1482
1483/********* Methode *********/
1484//++
1485void 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//++
1656void 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//++
1668void Histo::WriteSelf(POutPersist& s) const
1669//
1670// Ecriture ppersist
1671//--
1672{
1673char strg[256];
1674
1675// Que faut-il ecrire?
1676int errok = (err2) ? 1 : 0;
1677
1678// Ecriture entete pour identifier facilement
1679sprintf(strg,"bins=%6d NEnt=%15d errok=%1d",bins,nEntries,errok);
1680s.PutLine(strg);
1681sprintf(strg,"binw=%g min=%g max=%g",binWidth,min,max);
1682s.PutLine(strg);
1683sprintf(strg, "under=%g over=%g nHist=%g",under,over,nHist);
1684s.PutLine(strg);
1685
1686// Ecriture des valeurs
1687s.PutI4(bins);
1688s.PutI4(nEntries);
1689s.PutI4(errok);
1690
1691s.PutR4(binWidth);
1692s.PutR4(min);
1693s.PutR4(max);
1694s.PutR4(under);
1695s.PutR4(over);
1696
1697s.PutR8(nHist);
1698
1699// Ecriture des donnees Histo 1D
1700sprintf(strg,"Histo: Tableau des donnees %d",bins);
1701s.PutLine(strg);
1702s.PutR4s(data, bins);
1703
1704// Ecriture des erreurs
1705if(errok) {
1706 sprintf(strg,"Histo: Tableau des erreurs %d",bins);
1707 s.PutLine(strg);
1708 s.PutR8s(err2, bins);
1709}
1710
1711return;
1712}
1713
1714/********* Methode *********/
1715//++
1716void Histo::ReadSelf(PInPersist& s)
1717//
1718// Lecture ppersist
1719//--
1720{
1721char strg[256];
1722
1723Delete();
1724
1725// Lecture entete
1726s.GetLine(strg, 255);
1727s.GetLine(strg, 255);
1728s.GetLine(strg, 255);
1729
1730// Lecture des valeurs
1731s.GetI4(bins);
1732s.GetI4(nEntries);
1733int_4 errok;
1734s.GetI4(errok);
1735
1736s.GetR4(binWidth);
1737s.GetR4(min);
1738s.GetR4(max);
1739s.GetR4(under);
1740s.GetR4(over);
1741
1742s.GetR8(nHist);
1743
1744// Lecture des donnees Histo 1D
1745data = new float[bins];
1746s.GetLine(strg, 255);
1747s.GetR4s(data, bins);
1748
1749// Lecture des erreurs
1750if(errok) {
1751 s.GetLine(strg, 255);
1752 err2 = new double[bins];
1753 s.GetR8s(err2, bins);
1754}
1755
1756return;
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//--
Note: See TracBrowser for help on using the repository browser.