[220] | 1 | //
|
---|
[244] | 2 | // $Id: histos.cc,v 1.2 1999-04-22 16:18:36 ansari Exp $
|
---|
[220] | 3 | //
|
---|
| 4 |
|
---|
[244] | 5 | #include "machdefs.h"
|
---|
[220] | 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 | //--
|
---|