Changeset 1092 in Sophya for trunk/SophyaLib/HiStats/histos.cc
- Timestamp:
- Jul 26, 2000, 3:15:52 PM (25 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaLib/HiStats/histos.cc
r1064 r1092 1 1 // 2 // $Id: histos.cc,v 1. 9 2000-07-11 18:57:25 ansari Exp $2 // $Id: histos.cc,v 1.10 2000-07-26 13:15:15 ansari Exp $ 3 3 // 4 4 … … 22 22 /*! Constructeur par defaut */ 23 23 Histo::Histo() 24 : data(NULL), err2(NULL),25 under(0), over(0), nHist(0), nEntries(0),26 bins(0), min(0), max(0),24 : mData(NULL), mErr2(NULL), 25 mUnder(0), mOver(0), nHist(0), nEntries(0), 26 mBins(0), mMin(0), mMax(0), 27 27 binWidth(0) 28 28 { … … 32 32 /********* Methode *********/ 33 33 /*! Constructeur d'un histo de nBin bins allant de xMin a xMax */ 34 Histo::Histo(float xMin, float xMax, int nBin) 35 : data(new float[nBin]), err2(NULL), 36 under(0), over(0), nHist(0), nEntries(0), 37 bins(nBin), min(xMin), max(xMax), 38 binWidth((max - min)/nBin) 34 Histo::Histo(r_8 xMin, r_8 xMax, int_4 nBin) 35 : mData((nBin>0) ? new r_8[nBin] : NULL), 36 mErr2(NULL), 37 mUnder(0), mOver(0), nHist(0), nEntries(0), 38 mBins(nBin), mMin(xMin), mMax(xMax), 39 binWidth((nBin>0) ? (mMax-mMin)/nBin : 0) 39 40 { 40 41 Zero(); … … 43 44 44 45 /********* Methode *********/ 46 /*! Constructeur d'un histo de nBin bins allant de xMin a xMax */ 47 Histo::Histo(r_4 xMin, r_4 xMax, int_4 nBin) 48 : mData((nBin>0) ? new r_8[nBin] : NULL), 49 mErr2(NULL), 50 mUnder(0), mOver(0), nHist(0), nEntries(0), 51 mBins(nBin), mMin((r_8)xMin), mMax((r_8)xMax), 52 binWidth((nBin>0) ? (mMax-mMin)/nBin : 0) 53 { 54 Zero(); 55 END_CONSTRUCTOR 56 } 57 58 /********* Methode *********/ 45 59 /*! Constructeur par copie */ 46 60 Histo::Histo(const Histo& H) 47 : data((H.bins>0)? new float[H.bins] : NULL),48 err2(NULL),49 under(H.under), over(H.over), nHist(H.nHist), nEntries(H.nEntries),50 bins(H.bins), min(H.min), max(H.max),61 : mData((H.mBins>0)? new r_8[H.mBins] : NULL), 62 mErr2(NULL), 63 mUnder(H.mUnder), mOver(H.mOver), nHist(H.nHist), nEntries(H.nEntries), 64 mBins(H.mBins), mMin(H.mMin), mMax(H.mMax), 51 65 binWidth(H.binWidth) 52 66 { 53 if( bins>0) {54 memcpy( data, H.data, bins*sizeof(float));55 if( H. err2 != NULL ) {56 err2 = new double[bins];57 memcpy( err2, H.err2, bins*sizeof(double));67 if(mBins>0) { 68 memcpy(mData, H.mData, mBins*sizeof(r_8)); 69 if( H.mErr2 != NULL ) { 70 mErr2 = new r_8[mBins]; 71 memcpy(mErr2, H.mErr2, mBins*sizeof(r_8)); 58 72 } 59 73 } … … 65 79 void Histo::Delete() 66 80 { 67 if( data != NULL ) { delete[] data; data = NULL;}68 if( err2 != NULL ) { delete[] err2; err2 = NULL;}69 under = over = min = max = binWidth= 0.;81 if( mData != NULL ) { delete[] mData; mData = NULL;} 82 if( mErr2 != NULL ) { delete[] mErr2; mErr2 = NULL;} 83 mUnder = mOver = mMin = mMax = binWidth= 0.; 70 84 nHist = 0.; 71 nEntries = bins = 0;85 nEntries = mBins = 0; 72 86 } 73 87 … … 85 99 void Histo::Zero() 86 100 { 87 if( bins<=0) return;88 memset( data, 0, bins*sizeof(float));89 under = over = 0;101 if(mBins<=0) return; 102 memset(mData, 0, mBins*sizeof(r_8)); 103 mUnder = mOver = 0; 90 104 nHist = 0; 91 105 nEntries = 0; 92 if( err2 != NULL ) memset(err2, 0, bins*sizeof(double));106 if( mErr2 != NULL ) memset(mErr2, 0, mBins*sizeof(r_8)); 93 107 } 94 108 … … 99 113 void Histo::Errors() 100 114 { 101 if( bins<=0) return;102 if( err2==NULL) err2 = new double[bins];103 memset( err2, 0, bins*sizeof(double));115 if(mBins<=0) return; 116 if(mErr2==NULL) mErr2 = new r_8[mBins]; 117 memset(mErr2, 0, mBins*sizeof(r_8)); 104 118 } 105 119 … … 111 125 { 112 126 if(this == &h) return *this; 113 if( h. bins <= 0 ) {Delete(); return *this;}114 if( h. bins > bins ) Delete();115 if(!h. err2 && err2 ) { delete [] err2; err2=NULL;}116 if(! data) data = new float[h.bins];117 if(h. err2 && !err2 ) err2 = new double[h.bins];118 119 under = h.under;120 over = h.over;127 if( h.mBins <= 0 ) {Delete(); return *this;} 128 if( h.mBins > mBins ) Delete(); 129 if(!h.mErr2 && mErr2 ) { delete [] mErr2; mErr2=NULL;} 130 if(!mData) mData = new r_8[h.mBins]; 131 if(h.mErr2 && !mErr2 ) mErr2 = new r_8[h.mBins]; 132 133 mUnder = h.mUnder; 134 mOver = h.mOver; 121 135 nHist = h.nHist; 122 136 nEntries = h.nEntries; 123 bins = h.bins;124 m in = h.min;125 m ax = h.max;137 mBins = h.mBins; 138 mMin = h.mMin; 139 mMax = h.mMax; 126 140 binWidth = h.binWidth; 127 141 128 memcpy( data, h.data, bins*sizeof(float));129 if( err2) memcpy(err2, h.err2, bins*sizeof(double));142 memcpy(mData, h.mData, mBins*sizeof(r_8)); 143 if(mErr2) memcpy(mErr2, h.mErr2, mBins*sizeof(r_8)); 130 144 131 145 return *this; … … 136 150 Operateur de multiplication par une constante 137 151 */ 138 Histo& Histo::operator *= ( doubleb)139 { 140 doubleb2 = b*b;141 for(int i=0;i<bins;i++) {142 data[i] *= b;143 if( err2) err2[i] *= b2;144 } 145 under *= b;146 over *= b;152 Histo& Histo::operator *= (r_8 b) 153 { 154 r_8 b2 = b*b; 155 for(int_4 i=0;i<mBins;i++) { 156 mData[i] *= b; 157 if(mErr2) mErr2[i] *= b2; 158 } 159 mUnder *= b; 160 mOver *= b; 147 161 nHist *= b; 148 162 … … 153 167 Operateur de division par une constante 154 168 */ 155 Histo& Histo::operator /= ( doubleb)169 Histo& Histo::operator /= (r_8 b) 156 170 { 157 171 if (b==0.) THROW(inconsistentErr); 158 doubleb2 = b*b;159 for(int i=0;i<bins;i++) {160 data[i] /= b;161 if( err2) err2[i] /= b2;162 } 163 under /= b;164 over /= b;172 r_8 b2 = b*b; 173 for(int_4 i=0;i<mBins;i++) { 174 mData[i] /= b; 175 if(mErr2) mErr2[i] /= b2; 176 } 177 mUnder /= b; 178 mOver /= b; 165 179 nHist /= b; 166 180 … … 171 185 Operateur d'addition d'une constante 172 186 */ 173 Histo& Histo::operator += ( doubleb)174 { 175 for(int i=0;i<bins;i++) data[i] += b;176 under += b;177 over += b;178 nHist += bins*b;187 Histo& Histo::operator += (r_8 b) 188 { 189 for(int_4 i=0;i<mBins;i++) mData[i] += b; 190 mUnder += b; 191 mOver += b; 192 nHist += mBins*b; 179 193 180 194 return *this; … … 184 198 Operateur de soustraction d'une constante 185 199 */ 186 Histo& Histo::operator -= ( doubleb)187 { 188 for(int i=0;i<bins;i++) data[i] -= b;189 under -= b;190 over -= b;191 nHist -= bins*b;200 Histo& Histo::operator -= (r_8 b) 201 { 202 for(int_4 i=0;i<mBins;i++) mData[i] -= b; 203 mUnder -= b; 204 mOver -= b; 205 nHist -= mBins*b; 192 206 193 207 return *this; … … 200 214 Histo& Histo::operator += (const Histo& a) 201 215 { 202 if( bins!=a.bins) THROW(sizeMismatchErr);203 for(int i=0;i<bins;i++) {204 data[i] += a(i);205 if( err2 && a.err2) err2[i] += a.Error2(i);206 } 207 under += a.under;208 over += a.over;216 if(mBins!=a.mBins) THROW(sizeMismatchErr); 217 for(int_4 i=0;i<mBins;i++) { 218 mData[i] += a(i); 219 if(mErr2 && a.mErr2) mErr2[i] += a.Error2(i); 220 } 221 mUnder += a.mUnder; 222 mOver += a.mOver; 209 223 nHist += a.nHist; 210 224 nEntries += a.nEntries; … … 218 232 Histo& Histo::operator -= (const Histo& a) 219 233 { 220 if( bins!=a.bins) THROW(sizeMismatchErr);221 for(int i=0;i<bins;i++) {222 data[i] -= a(i);223 if( err2 && a.err2) err2[i] += a.Error2(i);224 } 225 under -= a.under;226 over -= a.over;234 if(mBins!=a.mBins) THROW(sizeMismatchErr); 235 for(int_4 i=0;i<mBins;i++) { 236 mData[i] -= a(i); 237 if(mErr2 && a.mErr2) mErr2[i] += a.Error2(i); 238 } 239 mUnder -= a.mUnder; 240 mOver -= a.mOver; 227 241 nHist -= a.nHist; 228 242 nEntries += a.nEntries; … … 236 250 Histo& Histo::operator *= (const Histo& a) 237 251 { 238 if( bins!=a.bins) THROW(sizeMismatchErr);252 if(mBins!=a.mBins) THROW(sizeMismatchErr); 239 253 nHist = 0.; 240 for(int i=0;i<bins;i++) {241 if( err2 && a.err2)242 err2[i] = a(i)*a(i)*err2[i] + data[i]*data[i]*a.Error2(i);243 data[i] *= a(i);244 nHist += data[i];245 } 246 under *= a.under;247 over *= a.over;254 for(int_4 i=0;i<mBins;i++) { 255 if(mErr2 && a.mErr2) 256 mErr2[i] = a(i)*a(i)*mErr2[i] + mData[i]*mData[i]*a.Error2(i); 257 mData[i] *= a(i); 258 nHist += mData[i]; 259 } 260 mUnder *= a.mUnder; 261 mOver *= a.mOver; 248 262 nEntries += a.nEntries; 249 263 … … 256 270 Histo& Histo::operator /= (const Histo& a) 257 271 { 258 if( bins!=a.bins) THROW(sizeMismatchErr);272 if(mBins!=a.mBins) THROW(sizeMismatchErr); 259 273 nHist = 0.; 260 for(int i=0;i<bins;i++) {274 for(int_4 i=0;i<mBins;i++) { 261 275 if(a(i)==0.) { 262 data[i]=0.;263 if( err2) err2[i]=0.;276 mData[i]=0.; 277 if(mErr2) mErr2[i]=0.; 264 278 continue; 265 279 } 266 if( err2 && a.err2)267 err2[i] = (err2[i] + data[i]/a(i)*data[i]/a(i)*a.Error2(i))280 if(mErr2 && a.mErr2) 281 mErr2[i] = (mErr2[i] + mData[i]/a(i)*mData[i]/a(i)*a.Error2(i)) 268 282 /(a(i)*a(i)); 269 data[i] /= a(i);270 nHist += data[i];271 } 272 if(a. under!=0.) under /= a.under; else under = 0.;273 if(a. over!=0.) over /= a.over; else over = 0.;283 mData[i] /= a(i); 284 nHist += mData[i]; 285 } 286 if(a.mUnder!=0.) mUnder /= a.mUnder; else mUnder = 0.; 287 if(a.mOver!=0.) mOver /= a.mOver; else mOver = 0.; 274 288 nEntries += a.nEntries; 275 289 … … 283 297 void Histo::GetAbsc(TVector<r_8> &v) 284 298 { 285 v.Realloc( bins);286 for(int i=0;i<bins;i++) v(i) = BinLowEdge(i);299 v.Realloc(mBins); 300 for(int_4 i=0;i<mBins;i++) v(i) = BinLowEdge(i); 287 301 return; 288 302 } … … 293 307 void Histo::GetValue(TVector<r_8> &v) 294 308 { 295 v.Realloc( bins);296 for(int i=0;i<bins;i++) v(i) = data[i];309 v.Realloc(mBins); 310 for(int_4 i=0;i<mBins;i++) v(i) = mData[i]; 297 311 return; 298 312 } … … 303 317 void Histo::GetError2(TVector<r_8> &v) 304 318 { 305 v.Realloc( bins);306 if(! err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}307 for(int i=0;i<bins;i++) v(i) = err2[i];319 v.Realloc(mBins); 320 if(!mErr2) {for(int_4 i=0;i<mBins;i++) v(i) = 0.; return;} 321 for(int_4 i=0;i<mBins;i++) v(i) = mErr2[i]; 308 322 return; 309 323 } … … 314 328 void Histo::GetError(TVector<r_8> &v) 315 329 { 316 v.Realloc( bins);317 if(! err2) {for(int i=0;i<bins;i++) v(i) = 0.; return;}318 for(int i=0;i<bins;i++) v(i) = Error(i);330 v.Realloc(mBins); 331 if(!mErr2) {for(int_4 i=0;i<mBins;i++) v(i) = 0.; return;} 332 for(int_4 i=0;i<mBins;i++) v(i) = Error(i); 319 333 return; 320 334 } … … 324 338 Remplissage du contenu de l'histo avec les valeurs d'un vecteur 325 339 */ 326 void Histo::PutValue(TVector<r_8> &v, int ierr)327 { 328 //if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);329 uint_4 n = (v.NElts()<(uint_4) bins) ? v.NElts(): (uint_4) bins;330 if(n>0) for( inti=0;i<n;i++) {331 data[i] = v(i);332 if( err2&&ierr) err2[i] = fabs(v(i));340 void Histo::PutValue(TVector<r_8> &v, int_4 ierr) 341 { 342 //if(v.NElts()<(uint_4) mBins) THROW(sizeMismatchErr); 343 uint_4 n = (v.NElts()<(uint_4) mBins) ? v.NElts(): (uint_4) mBins; 344 if(n>0) for(uint_4 i=0;i<n;i++) { 345 mData[i] = v(i); 346 if(mErr2&&ierr) mErr2[i] = fabs(v(i)); 333 347 } 334 348 return; … … 338 352 Addition du contenu de l'histo avec les valeurs d'un vecteur 339 353 */ 340 void Histo::PutValueAdd(TVector<r_8> &v, int ierr)341 { 342 //if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);343 uint_4 n = (v.NElts()<(uint_4) bins) ? v.NElts(): (uint_4) bins;344 if(n>0) for( inti=0;i<n;i++) {345 data[i] += v(i);346 if( err2 && ierr) err2[i] += fabs(v(i));354 void Histo::PutValueAdd(TVector<r_8> &v, int_4 ierr) 355 { 356 //if(v.NElts()<(uint_4) mBins) THROW(sizeMismatchErr); 357 uint_4 n = (v.NElts()<(uint_4) mBins) ? v.NElts(): (uint_4) mBins; 358 if(n>0) for(uint_4 i=0;i<n;i++) { 359 mData[i] += v(i); 360 if(mErr2 && ierr) mErr2[i] += fabs(v(i)); 347 361 } 348 362 return; … … 354 368 void Histo::PutError2(TVector<r_8> &v) 355 369 { 356 //if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);357 uint_4 n = (v.NElts()<(uint_4) bins) ? v.NElts(): (uint_4) bins;370 //if(v.NElts()<(uint_4) mBins) THROW(sizeMismatchErr); 371 uint_4 n = (v.NElts()<(uint_4) mBins) ? v.NElts(): (uint_4) mBins; 358 372 if(n>0) { 359 if(! err2) Errors();360 for( int i=0;i<n;i++) err2[i] = v(i);373 if(!mErr2) Errors(); 374 for(uint_4 i=0;i<n;i++) mErr2[i] = v(i); 361 375 } 362 376 return; … … 368 382 void Histo::PutError2Add(TVector<r_8> &v) 369 383 { 370 //if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);371 uint_4 n = (v.NElts()<(uint_4) bins) ? v.NElts(): (uint_4) bins;384 //if(v.NElts()<(uint_4) mBins) THROW(sizeMismatchErr); 385 uint_4 n = (v.NElts()<(uint_4) mBins) ? v.NElts(): (uint_4) mBins; 372 386 if(n>0) { 373 if(! err2) Errors();374 for( int i=0;i<n;i++) if(v(i)>0.) err2[i] += v(i);387 if(!mErr2) Errors(); 388 for(uint_4 i=0;i<n;i++) if(v(i)>0.) mErr2[i] += v(i); 375 389 } 376 390 return; … … 382 396 void Histo::PutError(TVector<r_8> &v) 383 397 { 384 //if(v.NElts()<(uint_4) bins) THROW(sizeMismatchErr);385 uint_4 n = (v.NElts()<(uint_4) bins) ? v.NElts(): (uint_4) bins;398 //if(v.NElts()<(uint_4) mBins) THROW(sizeMismatchErr); 399 uint_4 n = (v.NElts()<(uint_4) mBins) ? v.NElts(): (uint_4) mBins; 386 400 if(n>0) { 387 if(! err2) Errors();388 for( inti=0;i<n;i++)389 if(v(i)>0.) err2[i]=v(i)*v(i); else err2[i]=-v(i)*v(i);401 if(!mErr2) Errors(); 402 for(uint_4 i=0;i<n;i++) 403 if(v(i)>0.) mErr2[i]=v(i)*v(i); else mErr2[i]=-v(i)*v(i); 390 404 } 391 405 return; … … 396 410 Addition du contenu de l'histo pour abscisse x poids w 397 411 */ 398 void Histo::Add( float x, floatw)399 { 400 int numBin = FindBin(x);401 if (numBin<0) under += w;402 else if (numBin>= bins) over += w;412 void Histo::Add(r_8 x, r_8 w) 413 { 414 int_4 numBin = FindBin(x); 415 if (numBin<0) mUnder += w; 416 else if (numBin>=mBins) mOver += w; 403 417 else { 404 data[numBin] += w;405 if( err2!=NULL) err2[numBin] += w*w;418 mData[numBin] += w; 419 if(mErr2!=NULL) mErr2[numBin] += w*w; 406 420 nHist += w; 407 421 nEntries++; … … 413 427 Addition du contenu de l'histo bin numBin poids w 414 428 */ 415 void Histo::AddBin(int numBin, floatw)416 { 417 if (numBin<0) under += w;418 else if (numBin>= bins) over += w;429 void Histo::AddBin(int_4 numBin, r_8 w) 430 { 431 if (numBin<0) mUnder += w; 432 else if (numBin>=mBins) mOver += w; 419 433 else { 420 data[numBin] += w;421 if( err2!=NULL) err2[numBin] += w*w;434 mData[numBin] += w; 435 if(mErr2!=NULL) mErr2[numBin] += w*w; 422 436 nHist += w; 423 437 nEntries++; … … 429 443 Remplissage du contenu de l'histo pour abscisse x poids w 430 444 */ 431 void Histo::SetBin( float x, floatw)432 { 433 int numBin = FindBin(x);445 void Histo::SetBin(r_8 x, r_8 w) 446 { 447 int_4 numBin = FindBin(x); 434 448 SetBin(numBin,w); 435 449 } … … 439 453 Remplissage du contenu de l'histo pour numBin poids w 440 454 */ 441 void Histo::SetBin(int numBin, floatw)442 { 443 if (numBin<0) under = w;444 else if (numBin>= bins) over = w;455 void Histo::SetBin(int_4 numBin, r_8 w) 456 { 457 if (numBin<0) mUnder = w; 458 else if (numBin>=mBins) mOver = w; 445 459 else { 446 nHist -= data[numBin];447 data[numBin] = w;460 nHist -= mData[numBin]; 461 mData[numBin] = w; 448 462 nHist += w; 449 463 } … … 454 468 Remplissage des erreurs au carre pour abscisse x 455 469 */ 456 void Histo::SetErr2( float x, doublee2)457 { 458 int numBin = FindBin(x);470 void Histo::SetErr2(r_8 x, r_8 e2) 471 { 472 int_4 numBin = FindBin(x); 459 473 SetErr2(numBin,e2); 460 474 } … … 464 478 Remplissage des erreurs au carre pour numBin poids 465 479 */ 466 void Histo::SetErr2(int numBin, doublee2)467 { 468 if( err2==NULL) return;469 if ( numBin<0 || numBin>= bins ) return;470 err2[numBin] = e2;480 void Histo::SetErr2(int_4 numBin, r_8 e2) 481 { 482 if( mErr2==NULL) return; 483 if ( numBin<0 || numBin>=mBins ) return; 484 mErr2[numBin] = e2; 471 485 } 472 486 … … 475 489 Remplissage des erreurs pour abscisse x 476 490 */ 477 void Histo::SetErr( float x, floate)478 { 479 SetErr2(x, (double)e*e);491 void Histo::SetErr(r_8 x, r_8 e) 492 { 493 SetErr2(x, e*e); 480 494 } 481 495 … … 484 498 Remplissage des erreurs pour numBin 485 499 */ 486 void Histo::SetErr(int numBin, floate)487 { 488 SetErr2(numBin, (double)e*e);500 void Histo::SetErr(int_4 numBin, r_8 e) 501 { 502 SetErr2(numBin, e*e); 489 503 } 490 504 … … 493 507 Numero du bin ayant le contenu maximum 494 508 */ 495 int Histo::IMax() const496 { 497 int imx=0;498 if( bins==1) return imx;499 float mx=data[0];500 for (int i=1; i<bins; i++)501 if ( data[i]>mx) {imx = i; mx=data[i];}509 int_4 Histo::IMax() const 510 { 511 int_4 imx=0; 512 if(mBins==1) return imx; 513 r_8 mx=mData[0]; 514 for (int_4 i=1; i<mBins; i++) 515 if (mData[i]>mx) {imx = i; mx=mData[i];} 502 516 return imx; 503 517 } … … 507 521 Numero du bin ayant le contenu minimum 508 522 */ 509 int Histo::IMin() const510 { 511 int imx=0;512 if( bins==1) return imx;513 float mx=data[0];514 for (int i=1; i<bins; i++)515 if ( data[i]<mx) {imx = i; mx=data[i];}523 int_4 Histo::IMin() const 524 { 525 int_4 imx=0; 526 if(mBins==1) return imx; 527 r_8 mx=mData[0]; 528 for (int_4 i=1; i<mBins; i++) 529 if (mData[i]<mx) {imx = i; mx=mData[i];} 516 530 return imx; 517 531 } … … 521 535 Valeur le contenu maximum 522 536 */ 523 floatHisto::VMax() const524 { 525 float mx=data[0];526 if( bins==1) return mx;527 for (int i=1;i<bins;i++) if(data[i]>mx) mx=data[i];537 r_8 Histo::VMax() const 538 { 539 r_8 mx=mData[0]; 540 if(mBins==1) return mx; 541 for (int_4 i=1;i<mBins;i++) if(mData[i]>mx) mx=mData[i]; 528 542 return mx; 529 543 } … … 533 547 Valeur le contenu minimum 534 548 */ 535 floatHisto::VMin() const536 { 537 float mx=data[0];538 if( bins==1) return mx;539 for (int i=1;i<bins;i++) if(data[i]<mx) mx=data[i];549 r_8 Histo::VMin() const 550 { 551 r_8 mx=mData[0]; 552 if(mBins==1) return mx; 553 for (int_4 i=1;i<mBins;i++) if(mData[i]<mx) mx=mData[i]; 540 554 return mx; 541 555 } … … 545 559 Valeur moyenne 546 560 */ 547 floatHisto::Mean() const548 { 549 doublen = 0;550 doublesx = 0;551 for (int i=0; i<bins; i++) {552 double v = (data[i]>=0.) ? data[i] : -data[i];561 r_8 Histo::Mean() const 562 { 563 r_8 n = 0; 564 r_8 sx = 0; 565 for (int_4 i=0; i<mBins; i++) { 566 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 553 567 n += v; 554 568 sx += BinCenter(i)*v; 555 569 } 556 if(n>0.) return sx/n; else return m in;570 if(n>0.) return sx/n; else return mMin; 557 571 } 558 572 … … 561 575 Valeur du sigma 562 576 */ 563 floatHisto::Sigma() const564 { 565 doublen = 0;566 doublesx = 0;567 doublesx2= 0;568 for (int i=0; i<bins; i++) {569 double v = (data[i]>=0.) ? data[i] : -data[i];577 r_8 Histo::Sigma() const 578 { 579 r_8 n = 0; 580 r_8 sx = 0; 581 r_8 sx2= 0; 582 for (int_4 i=0; i<mBins; i++) { 583 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 570 584 n += v; 571 585 sx += BinCenter(i)*v; … … 583 597 Valeur de la moyenne calculee entre les bin il et ih 584 598 */ 585 float Histo::MeanLH(int il,intih) const586 { 587 if( ih < il ) { il = 0; ih = bins-1; }599 r_8 Histo::MeanLH(int_4 il,int_4 ih) const 600 { 601 if( ih < il ) { il = 0; ih = mBins-1; } 588 602 if( il < 0 ) il = 0; 589 if( ih >= bins ) ih = bins-1;590 doublen = 0;591 doublesx = 0;592 for (int i=il; i<=ih; i++) {593 double v = (data[i]>=0.) ? data[i] : -data[i];603 if( ih >= mBins ) ih = mBins-1; 604 r_8 n = 0; 605 r_8 sx = 0; 606 for (int_4 i=il; i<=ih; i++) { 607 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 594 608 n += v; 595 609 sx += BinCenter(i)*v; … … 602 616 Valeur de la moyenne calculee entre les bin il et ih 603 617 */ 604 float Histo::SigmaLH(int il,intih) const605 { 606 if( ih < il ) { il = 0; ih = bins - 1; }618 r_8 Histo::SigmaLH(int_4 il,int_4 ih) const 619 { 620 if( ih < il ) { il = 0; ih = mBins - 1; } 607 621 if( il < 0 ) il = 0; 608 if( ih >= bins ) ih = bins - 1;609 doublen = 0;610 doublesx = 0;611 doublesx2= 0;612 for (int i=il; i<=ih; i++) {613 double v = (data[i]>=0.) ? data[i] : -data[i];622 if( ih >= mBins ) ih = mBins - 1; 623 r_8 n = 0; 624 r_8 sx = 0; 625 r_8 sx2= 0; 626 for (int_4 i=il; i<=ih; i++) { 627 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 614 628 n += v; 615 629 sx += BinCenter(i)*v; … … 626 640 Valeur de la moyenne calculee entre x0-dx et x0+dx 627 641 */ 628 float Histo::Mean(float x0, floatdx) const629 { 630 doublesdata = 0;631 doublesx = 0;632 int iMin = FindBin(x0-dx);642 r_8 Histo::Mean(r_8 x0, r_8 dx) const 643 { 644 r_8 sdata = 0; 645 r_8 sx = 0; 646 int_4 iMin = FindBin(x0-dx); 633 647 if (iMin<0) iMin = 0; 634 int iMax = FindBin(x0+dx);635 if (iMax> bins) iMax=bins;636 for (int i=iMin; i<iMax; i++) {637 double v = (data[i]>=0.) ? data[i] : -data[i];648 int_4 iMax = FindBin(x0+dx); 649 if (iMax>mBins) iMax=mBins; 650 for (int_4 i=iMin; i<iMax; i++) { 651 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 638 652 sx += BinCenter(i)*v; 639 653 sdata += v; 640 654 } 641 if(sdata>0.) return sx/sdata; else return m in;655 if(sdata>0.) return sx/sdata; else return mMin; 642 656 } 643 657 … … 646 660 Valeur du sigma calcule entre x0-dx et x0+dx 647 661 */ 648 float Histo::Sigma(float x0, floatdx) const649 { 650 doublesx = 0;651 doublesx2= 0;652 doublesdata = 0;653 int iMin = FindBin(x0-dx);662 r_8 Histo::Sigma(r_8 x0, r_8 dx) const 663 { 664 r_8 sx = 0; 665 r_8 sx2= 0; 666 r_8 sdata = 0; 667 int_4 iMin = FindBin(x0-dx); 654 668 if (iMin<0) iMin = 0; 655 int iMax = FindBin(x0+dx);656 if (iMax> bins) iMax=bins;657 for (int i=iMin; i<iMax; i++) {658 double v = (data[i]>=0.) ? data[i] : -data[i];669 int_4 iMax = FindBin(x0+dx); 670 if (iMax>mBins) iMax=mBins; 671 for (int_4 i=iMin; i<iMax; i++) { 672 r_8 v = (mData[i]>=0.) ? mData[i] : -mData[i]; 659 673 sx += BinCenter(i)*v; 660 674 sx2+= BinCenter(i)*BinCenter(i)*v; … … 669 683 Valeur de la moyenne et du sigma nettoyes 670 684 */ 671 float Histo::CleanedMean(float& sigma) const685 r_8 Histo::CleanedMean(r_8& sigma) const 672 686 { 673 687 if (!nHist) return 0; 674 688 // on fait ca de facon bete, on pourra raffiner apres... 675 floatx0 = Mean();676 floats = Sigma()+binWidth;689 r_8 x0 = Mean(); 690 r_8 s = Sigma()+binWidth; 677 691 678 for (int i=0; i<3; i++) {679 floatxx0 = Mean(x0,5*s);692 for (int_4 i=0; i<3; i++) { 693 r_8 xx0 = Mean(x0,5*s); 680 694 s = Sigma(x0,5*s)+binWidth; 681 695 x0 = xx0; … … 689 703 Valeur de la moyenne nettoyee 690 704 */ 691 floatHisto::CleanedMean() const705 r_8 Histo::CleanedMean() const 692 706 { 693 707 if (!nHist) return 0; 694 floats=0;708 r_8 s=0; 695 709 return CleanedMean(s); 696 710 } … … 700 714 Retourne le nombre de bins non-nul 701 715 */ 702 int Histo::BinNonNul() const703 { 704 int non=0;705 for (int i=0;i<bins;i++) if( data[i] != 0. ) non++;716 int_4 Histo::BinNonNul() const 717 { 718 int_4 non=0; 719 for (int_4 i=0;i<mBins;i++) if( mData[i] != 0. ) non++; 706 720 return non; 707 721 } … … 711 725 Retourne le nombre de bins ayant une erreur non-nulle 712 726 */ 713 int Histo::ErrNonNul() const714 { 715 if( err2==NULL) return -1;716 int non=0;717 for (int i=0;i<bins;i++) if( err2[i] != 0. ) non++;727 int_4 Histo::ErrNonNul() const 728 { 729 if(mErr2==NULL) return -1; 730 int_4 non=0; 731 for (int_4 i=0;i<mBins;i++) if( mErr2[i] != 0. ) non++; 718 732 return non; 719 733 } … … 724 738 entre le bin 0 et ce bin (y compris le contenu de ce bin). 725 739 */ 726 int Histo::BinPercent(floatper) const727 { 728 doublen = nHist*per;729 doubles = 0.;730 int i;731 for(i=0; i< bins; i++ ) {732 s += data[i];740 int_4 Histo::BinPercent(r_8 per) const 741 { 742 r_8 n = nHist*per; 743 r_8 s = 0.; 744 int_4 i; 745 for(i=0; i<mBins; i++ ) { 746 s += mData[i]; 733 747 if( s >= n ) break; 734 748 } 735 if( i == bins ) i = bins-1;749 if( i == mBins ) i = mBins-1; 736 750 return i; 737 751 } … … 747 761 imax = bin tel que nombre d'entrees entre I et imax = N2 * per 748 762 Return: <0 si echec 749 m in(I-imin,imax-I) si ok763 mMin(I-imin,imax-I) si ok 750 764 \endverbatim 751 765 */ 752 int Histo::BinPercent(float x,floatper,int& imin,int& imax)766 int_4 Histo::BinPercent(r_8 x,r_8 per,int& imin,int& imax) 753 767 { 754 768 imin = imax = -1; 755 769 if( per <= 0. ) return -1; 756 770 757 int I = FindBin(x);758 if( I<0 || I>= bins ) return -2;759 760 double N1 = 0.; for(int i=0; i<=I; i++) N1 += data[i]; N1 *= per;761 double N2 = 0.; {for(int i=I; i<bins; i++) N2 += data[i]; N2 *= per;}762 763 doublen = 0.;764 for(imin=I; imin>=0; imin--) { n += data[imin]; if(n>=N1) break; }771 int_4 I = FindBin(x); 772 if( I<0 || I>=mBins ) return -2; 773 774 r_8 N1 = 0.; for(int_4 i=0; i<=I; i++) N1 += mData[i]; N1 *= per; 775 r_8 N2 = 0.; {for(int_4 i=I; i<mBins; i++) N2 += mData[i]; N2 *= per;} 776 777 r_8 n = 0.; 778 for(imin=I; imin>=0; imin--) { n += mData[imin]; if(n>=N1) break; } 765 779 if( imin<0 ) imin = 0; 766 780 // cout<<"I="<<I<<" imin="<<imin<<" n="<<n<<" N1="<<N1<<endl; 767 781 768 782 n = 0.; 769 for(imax=I; imax< bins; imax++) { n += data[imax]; if(n>=N2) break; }770 if( imax>= bins ) imax = bins-1;783 for(imax=I; imax<mBins; imax++) { n += mData[imax]; if(n>=N2) break; } 784 if( imax>=mBins ) imax = mBins-1; 771 785 // cout<<"I="<<I<<" imax="<<imax<<" n="<<n<<" N2="<<N2<<endl; 772 786 … … 778 792 Idem precedent mais renvoie xmin et xmax 779 793 */ 780 int Histo::BinPercent(float x,float per,float& xmin,float& xmax)794 int_4 Histo::BinPercent(r_8 x,r_8 per,r_8& xmin,r_8& xmax) 781 795 { 782 796 xmin = xmax = 0.; 783 int imin,imax;784 int rc = BinPercent(x,per,imin,imax);797 int_4 imin,imax; 798 int_4 rc = BinPercent(x,per,imin,imax); 785 799 if( rc >= 0 ) { xmin = BinLowEdge(imin); xmax = BinHighEdge(imax); } 786 800 return rc; … … 792 806 si norm <= 0 : pas de normalisation, integration seule 793 807 */ 794 void Histo::HInteg( floatnorm)795 { 796 if( bins <= 0 ) return; // createur par default797 if( bins>1)798 for(int i=1; i<bins; i++) {799 data[i] += data[i-1];800 if( err2!=NULL) err2[i] += err2[i-1];808 void Histo::HInteg(r_8 norm) 809 { 810 if( mBins <= 0 ) return; // createur par default 811 if(mBins>1) 812 for(int_4 i=1; i<mBins; i++) { 813 mData[i] += mData[i-1]; 814 if(mErr2!=NULL) mErr2[i] += mErr2[i-1]; 801 815 } 802 816 if( norm <= 0. ) return; 803 norm /= data[bins-1];804 for(int i=0; i<bins; i++) {805 data[i] *= norm;806 if( err2!=NULL) err2[i] *= norm*norm;817 norm /= mData[mBins-1]; 818 for(int_4 i=0; i<mBins; i++) { 819 mData[i] *= norm; 820 if(mErr2!=NULL) mErr2[i] *= norm*norm; 807 821 } 808 822 } … … 814 828 void Histo::HDeriv() 815 829 { 816 if( bins <= 0 ) return; // createur par default817 if( bins <= 1 ) return;818 float *temp = new float[bins];819 memcpy(temp, data, bins*sizeof(float));820 if( bins>=3) for(int i=1; i<bins-1; i++)821 temp[i] = ( data[i+1]-data[i-1])/(2.*binWidth);822 temp[0] = ( data[1]-data[0])/binWidth;823 temp[ bins-1] = (data[bins-1]-data[bins-2])/binWidth;824 memcpy( data, temp, bins*sizeof(float));830 if( mBins <= 0 ) return; // createur par default 831 if( mBins <= 1 ) return; 832 r_8 *temp = new r_8[mBins]; 833 memcpy(temp, mData, mBins*sizeof(r_8)); 834 if(mBins>=3) for(int_4 i=1; i<mBins-1; i++) 835 temp[i] = (mData[i+1]-mData[i-1])/(2.*binWidth); 836 temp[0] = (mData[1]-mData[0])/binWidth; 837 temp[mBins-1] = (mData[mBins-1]-mData[mBins-2])/binWidth; 838 memcpy(mData, temp, mBins*sizeof(r_8)); 825 839 delete [] temp; 826 840 } … … 830 844 Pour rebinner l'histogramme sur nbinew bins 831 845 */ 832 void Histo::HRebin(int nbinew)833 { 834 if( bins <= 0 ) return; // createur par default846 void Histo::HRebin(int_4 nbinew) 847 { 848 if( mBins <= 0 ) return; // createur par default 835 849 if( nbinew <= 0 ) return; 836 850 … … 839 853 840 854 // Le nombre de bins est il plus grand ?? 841 if( nbinew > bins ) {842 delete [] data; data = NULL;843 data = new float[nbinew];844 if( err2 != NULL ) {845 delete [] err2; err2 = NULL;846 err2 = new double[nbinew];855 if( nbinew > mBins ) { 856 delete [] mData; mData = NULL; 857 mData = new r_8[nbinew]; 858 if( mErr2 != NULL ) { 859 delete [] mErr2; mErr2 = NULL; 860 mErr2 = new r_8[nbinew]; 847 861 } 848 862 } 849 863 850 864 // remise en forme de this 851 bins = nbinew;865 mBins = nbinew; 852 866 this->Zero(); 853 binWidth = (m ax-min)/bins;867 binWidth = (mMax-mMin)/mBins; 854 868 // On recopie les parametres qui n'ont pas changes 855 under = H.under;856 over = H.over;869 mUnder = H.mUnder; 870 mOver = H.mOver; 857 871 858 872 859 873 // Remplissage 860 for(int i=0;i<bins;i++) {861 floatxmi = BinLowEdge(i);862 floatxma = BinHighEdge(i);863 int imi = H.FindBin(xmi);874 for(int_4 i=0;i<mBins;i++) { 875 r_8 xmi = BinLowEdge(i); 876 r_8 xma = BinHighEdge(i); 877 int_4 imi = H.FindBin(xmi); 864 878 if( imi < 0 ) imi = 0; 865 int ima = H.FindBin(xma);866 if( ima >= H. bins ) ima = H.bins-1;867 doublew = 0.;879 int_4 ima = H.FindBin(xma); 880 if( ima >= H.mBins ) ima = H.mBins-1; 881 r_8 w = 0.; 868 882 if( ima == imi ) { 869 883 w = H(imi) * binWidth/H.binWidth; … … 872 886 w += H(ima) * (xma -H.BinLowEdge(ima))/H.binWidth; 873 887 if( ima > imi+1 ) 874 for(int ii=imi+1;ii<ima;ii++) w += H(ii);888 for(int_4 ii=imi+1;ii<ima;ii++) w += H(ii); 875 889 } 876 AddBin(i, (float)w);890 AddBin(i, w); 877 891 } 878 892 … … 887 901 consecutifs, le bin le plus a droite est pris. 888 902 */ 889 int Histo::MaxiLocal(float& maxi,int& imax,float& maxn,int& imaxn)890 { 891 int nml = 0;903 int_4 Histo::MaxiLocal(r_8& maxi,int& imax,r_8& maxn,int& imaxn) 904 { 905 int_4 nml = 0; 892 906 imax = imaxn = -1; 893 maxi = maxn = -1. f;907 maxi = maxn = -1.; 894 908 895 909 bool up = true; 896 910 bool down = false; 897 for(int i=0;i<bins;i++) {898 if( !up ) if( data[i] > data[i-1] ) up = true;911 for(int_4 i=0;i<mBins;i++) { 912 if( !up ) if( mData[i] > mData[i-1] ) up = true; 899 913 if( up && !down ) { 900 if( i == bins-1 ) down=true;901 else if( data[i] > data[i+1] ) down=true;914 if( i == mBins-1 ) down=true; 915 else if( mData[i] > mData[i+1] ) down=true; 902 916 } 903 917 … … 906 920 up = down = false; 907 921 if( imax < 0 ) { 908 imax = i; maxi = data[i];909 } else if( data[i] >= maxi ) {922 imax = i; maxi = mData[i]; 923 } else if( mData[i] >= maxi ) { 910 924 imaxn = imax; maxn = maxi; 911 imax = i; maxi = data[i];925 imax = i; maxi = mData[i]; 912 926 } else { 913 if( imaxn < 0 || data[i] >= maxn ) { imaxn = i; maxn = data[i]; }927 if( imaxn < 0 || mData[i] >= maxn ) { imaxn = i; maxn = mData[i]; } 914 928 } 915 929 } … … 925 939 L'histo est suppose etre remplit de valeurs positives 926 940 */ 927 float Histo::FitMax(int degree, float frac, intdebug) const941 r_8 Histo::FitMax(int_4 degree, r_8 frac, int_4 debug) const 928 942 { 929 943 if (degree < 2 || degree > 3) degree = 3; … … 935 949 if (NEntries() < 1) THROW(inconsistentErr); 936 950 937 int iMax = IMax();938 floathmax = (*this)(iMax);939 floatxCenter = BinCenter(iMax);951 int_4 iMax = IMax(); 952 r_8 hmax = (*this)(iMax); 953 r_8 xCenter = BinCenter(iMax); 940 954 941 955 if(debug>1) … … 944 958 /* find limits at frac*hmax */ 945 959 946 floatlimit = frac*hmax;947 948 volatile int iLow = iMax;960 r_8 limit = frac*hmax; 961 962 volatile int_4 iLow = iMax; 949 963 while (iLow>0 && (*this)(iLow)>limit) iLow--; 950 964 951 volatile int iHigh = iMax;965 volatile int_4 iHigh = iMax; 952 966 while (iHigh<NBins()-1 && (*this)(iHigh)>limit) iHigh++; 953 967 954 int nLowHigh;968 int_4 nLowHigh; 955 969 for(;;) { 956 970 nLowHigh = 0; 957 for (int i=iLow; i<=iHigh; i++) if((*this)(i)>0) {958 if(! err2) nLowHigh++;971 for (int_4 i=iLow; i<=iHigh; i++) if((*this)(i)>0) { 972 if(!mErr2) nLowHigh++; 959 973 else if(Error2(i)>0.) nLowHigh++; 960 974 } … … 977 991 TVector<r_8> e2Fit(nLowHigh); 978 992 TVector<r_8> errcoef(1); 979 int ii = 0;980 for (int j=iLow; j<=iHigh; j++) {993 int_4 ii = 0; 994 for (int_4 j=iLow; j<=iHigh; j++) { 981 995 if ((*this)(j)>0) { 982 if(! err2) {996 if(!mErr2) { 983 997 xFit(ii) = BinCenter(j)-xCenter; 984 998 yFit(ii) = (*this)(j); … … 994 1008 } 995 1009 if(debug>4) { 996 int k;1010 int_4 k; 997 1011 for(k=0;k<nLowHigh;k++) cout<<" "<<xFit(k); cout<<endl; 998 1012 for(k=0;k<nLowHigh;k++) cout<<" "<<yFit(k); cout<<endl; … … 1009 1023 } ENDTRY 1010 1024 1011 int DPolDeg = pol.Degre();1012 floatfd = 0.;1025 int_4 DPolDeg = pol.Degre(); 1026 r_8 fd = 0.; 1013 1027 1014 1028 if (DPolDeg == 0) { … … 1026 1040 } else if (DPolDeg == 1) { 1027 1041 // on est dans le cas d'un fit de parabole 1028 doubler=0;1042 r_8 r=0; 1029 1043 if (pol.Root1(r)==0) THROW(inconsistentErr); 1030 1044 fd = r + xCenter; 1031 1045 } else if (DPolDeg == 2) { 1032 1046 // on est dans le cas d'un fit de cubique 1033 doubler1=0;1034 doubler2=0;1047 r_8 r1=0; 1048 r_8 r2=0; 1035 1049 if (pol.Root2(r1,r2) == 0) THROW(inconsistentErr); 1036 1050 pol.Derivate(); … … 1041 1055 } 1042 1056 1043 if(fd>m ax) fd = max;1044 if(fd<m in) fd = min;1057 if(fd>mMax) fd = mMax; 1058 if(fd<mMin) fd = mMin; 1045 1059 1046 1060 if (debug>1) … … 1060 1074 L'histo est suppose etre remplit de valeurs positives 1061 1075 */ 1062 float Histo::FindWidth(float frac, intdebug) const1063 { 1064 floatxmax = BinCenter( IMax() );1076 r_8 Histo::FindWidth(r_8 frac, int_4 debug) const 1077 { 1078 r_8 xmax = BinCenter( IMax() ); 1065 1079 return FindWidth(xmax,frac,debug); 1066 1080 } … … 1071 1085 L'histo est suppose etre remplit de valeurs positives 1072 1086 */ 1073 float Histo::FindWidth(float xmax,float frac, intdebug) const1087 r_8 Histo::FindWidth(r_8 xmax,r_8 frac, int_4 debug) const 1074 1088 { 1075 1089 if (frac <= 0 || frac >= 1.) frac = 0.5; … … 1085 1099 if (NBins() < 3) THROW(inconsistentErr); 1086 1100 1087 int iMax = FindBin(xmax);1101 int_4 iMax = FindBin(xmax); 1088 1102 if (iMax<0 || iMax>=NBins()) THROW(inconsistentErr); 1089 float hmax = data[iMax];1090 floatlimit = frac*hmax;1103 r_8 hmax = mData[iMax]; 1104 r_8 limit = frac*hmax; 1091 1105 if (debug > 1) 1092 1106 cout << " Max histo[" << iMax << "] = " << hmax 1093 1107 << ", limite " << limit << endl; 1094 1108 1095 int iLow = iMax;1096 while (iLow>=0 && data[iLow]>limit) iLow--;1109 int_4 iLow = iMax; 1110 while (iLow>=0 && mData[iLow]>limit) iLow--; 1097 1111 if( iLow < 0 ) iLow = 0; 1098 1112 1099 int iHigh = iMax;1100 while (iHigh<NBins() && data[iHigh]>limit) iHigh++;1113 int_4 iHigh = iMax; 1114 while (iHigh<NBins() && mData[iHigh]>limit) iHigh++; 1101 1115 if( iHigh >=NBins() ) iHigh = NBins()-1; 1102 1116 1103 floatxlow = BinCenter(iLow);1104 float ylow = data[iLow];1105 1106 floatxlow1=xlow, ylow1=ylow;1117 r_8 xlow = BinCenter(iLow); 1118 r_8 ylow = mData[iLow]; 1119 1120 r_8 xlow1=xlow, ylow1=ylow; 1107 1121 if(iLow+1<NBins()) { 1108 1122 xlow1 = BinCenter(iLow+1); 1109 ylow1 = data[iLow+1];1110 } 1111 1112 floatxhigh = BinCenter(iHigh);1113 float yhigh = data[iHigh];1114 1115 floatxhigh1=xhigh, yhigh1=yhigh;1123 ylow1 = mData[iLow+1]; 1124 } 1125 1126 r_8 xhigh = BinCenter(iHigh); 1127 r_8 yhigh = mData[iHigh]; 1128 1129 r_8 xhigh1=xhigh, yhigh1=yhigh; 1116 1130 if(iHigh-1>=0) { 1117 1131 xhigh1 = BinCenter(iHigh-1); 1118 yhigh1 = data[iHigh-1];1119 } 1120 1121 floatxlpred,xhpred,wd;1132 yhigh1 = mData[iHigh-1]; 1133 } 1134 1135 r_8 xlpred,xhpred,wd; 1122 1136 1123 1137 if(ylow1>ylow) { … … 1148 1162 Cf suivant mais im est le bin du maximum de l'histo 1149 1163 */ 1150 int Histo::EstimeMax(float& xm,intSzPav)1151 { 1152 int im = IMax();1164 int_4 Histo::EstimeMax(r_8& xm,int_4 SzPav) 1165 { 1166 int_4 im = IMax(); 1153 1167 return EstimeMax(im,xm,SzPav); 1154 1168 } … … 1166 1180 \endverbatim 1167 1181 */ 1168 int Histo::EstimeMax(int& im,float& xm,intSzPav)1182 int_4 Histo::EstimeMax(int& im,r_8& xm,int_4 SzPav) 1169 1183 { 1170 1184 xm = 0; 1171 1185 if( SzPav <= 0 ) return -1; 1172 if( im < 0 || im >= bins ) return -1;1186 if( im < 0 || im >= mBins ) return -1; 1173 1187 1174 1188 if( SzPav%2 == 0 ) SzPav++; 1175 1189 SzPav = (SzPav-1)/2; 1176 1190 1177 int rc = 0;1178 doubledxm = 0, wx = 0;1179 for(int i=im-SzPav;i<=im+SzPav;i++) {1180 if( i<0 || i>= bins ) {rc=1; continue;}1191 int_4 rc = 0; 1192 r_8 dxm = 0, wx = 0; 1193 for(int_4 i=im-SzPav;i<=im+SzPav;i++) { 1194 if( i<0 || i>= mBins ) {rc=1; continue;} 1181 1195 dxm += BinCenter(i) * (*this)(i); 1182 1196 wx += (*this)(i); … … 1198 1212 et a droite (widthD) 1199 1213 */ 1200 void Histo::EstimeWidthS( float frac,float& widthG,float& widthD)1201 { 1202 int i;1214 void Histo::EstimeWidthS(r_8 frac,r_8& widthG,r_8& widthD) 1215 { 1216 int_4 i; 1203 1217 widthG = widthD = -1.; 1204 if( bins<=1 || frac<=0. || frac>=1. ) return;1205 1206 int imax = 0;1207 float maxi = data[0];1208 for(i=1;i< bins;i++) if(data[i]>maxi) {imax=i; maxi=data[i];}1209 floatxmax = BinCenter(imax);1210 floatmaxf = maxi * frac;1218 if( mBins<=1 || frac<=0. || frac>=1. ) return; 1219 1220 int_4 imax = 0; 1221 r_8 maxi = mData[0]; 1222 for(i=1;i<mBins;i++) if(mData[i]>maxi) {imax=i; maxi=mData[i];} 1223 r_8 xmax = BinCenter(imax); 1224 r_8 maxf = maxi * frac; 1211 1225 1212 1226 // recherche du sigma a gauche 1213 1227 widthG = 0.; 1214 for(i=imax;i>=0;i--) if( data[i] <= maxf ) break;1228 for(i=imax;i>=0;i--) if( mData[i] <= maxf ) break; 1215 1229 if(i<0) i=0; 1216 1230 if(i<imax) { 1217 if( data[i+1] != data[i] ) {1231 if( mData[i+1] != mData[i] ) { 1218 1232 widthG = BinCenter(i) + binWidth 1219 * (maxf- data[i])/(data[i+1]-data[i]);1233 * (maxf-mData[i])/(mData[i+1]-mData[i]); 1220 1234 widthG = xmax - widthG; 1221 1235 } else widthG = xmax - BinCenter(i); … … 1224 1238 // recherche du sigma a droite 1225 1239 widthD = 0.; 1226 for(i=imax;i< bins;i++) if( data[i] <= maxf ) break;1227 if(i>= bins) i=bins-1;1240 for(i=imax;i<mBins;i++) if( mData[i] <= maxf ) break; 1241 if(i>=mBins) i=mBins-1; 1228 1242 if(i>imax) { 1229 if( data[i] != data[i-1] ) {1243 if( mData[i] != mData[i-1] ) { 1230 1244 widthD = BinCenter(i) - binWidth 1231 * (maxf- data[i])/(data[i-1]-data[i]);1245 * (maxf-mData[i])/(mData[i-1]-mData[i]); 1232 1246 widthD -= xmax; 1233 1247 } else widthD = BinCenter(i) - xmax; … … 1257 1271 \endverbatim 1258 1272 */ 1259 int Histo::Fit(GeneralFit& gfit,unsigned short typ_err)1273 int_4 Histo::Fit(GeneralFit& gfit,unsigned short typ_err) 1260 1274 { 1261 1275 if(NBins()<=0) return -1000; … … 1264 1278 GeneralFitData mydata(1,NBins()); 1265 1279 1266 for(int i=0;i<NBins();i++) {1267 double x = (double)BinCenter(i);1268 double f = (double)(*this)(i);1269 double saf = sqrt(fabs((double)f)); if(saf<1.) saf=1.;1270 doublee=0.;1280 for(int_4 i=0;i<NBins();i++) { 1281 r_8 x = BinCenter(i); 1282 r_8 f = (*this)(i); 1283 r_8 saf = sqrt(fabs( f)); if(saf<1.) saf=1.; 1284 r_8 e=0.; 1271 1285 if(typ_err==0) {if(HasErrors()) e=Error(i); else e=1.;} 1272 1286 else if(typ_err==1) {if(HasErrors()) e=Error(i); else e=saf;} … … 1295 1309 TVector<r_8> par = gfit.GetParm(); 1296 1310 Histo h(*this); 1297 for(int i=0;i<NBins();i++) {1298 double x = (double)BinCenter(i);1299 h(i) -= (float)f->Value(&x,par.Data());1311 for(int_4 i=0;i<NBins();i++) { 1312 r_8 x = BinCenter(i); 1313 h(i) -= f->Value(&x,par.Data()); 1300 1314 } 1301 1315 return h; … … 1314 1328 TVector<r_8> par = gfit.GetParm(); 1315 1329 Histo h(*this); 1316 for(int i=0;i<NBins();i++) {1317 double x = (double)BinCenter(i);1318 h(i) = (float)f->Value(&x,par.Data());1330 for(int_4 i=0;i<NBins();i++) { 1331 r_8 x = BinCenter(i); 1332 h(i) = f->Value(&x,par.Data()); 1319 1333 } 1320 1334 return h; … … 1334 1348 sinon : hmin hmax 1335 1349 pflag < 0 : on imprime les informations (nbin,min,...) sans l'histogramme 1336 = 0 : on imprime "BinCenter(i) data[i]" (note "... ...")1350 = 0 : on imprime "BinCenter(i) mData[i]" (note "... ...") 1337 1351 bit 0 on : (v=1): numero_bin "... ..." 1338 1352 bit 1 on : (v=2): "... ..." erreur 1339 1353 \endverbatim 1340 1354 */ 1341 void Histo::PrintF(FILE * fp, int hdyn,float hmin, float hmax,intpflag,1342 int il, intih)1355 void Histo::PrintF(FILE * fp, int_4 hdyn,r_8 hmin, r_8 hmax,int_4 pflag, 1356 int_4 il, int_4 ih) 1343 1357 { 1344 1358 1345 if( il > ih ) { il = 0; ih = bins-1; }1359 if( il > ih ) { il = 0; ih = mBins-1; } 1346 1360 if( il < 0 ) il = 0; 1347 if( ih >= bins ) ih = bins-1;1348 1349 double dhmin = (double)hmin;1350 double dhmax = (double)hmax;1351 doublehb,hbmn,hbmx;1361 if( ih >= mBins ) ih = mBins-1; 1362 1363 r_8 dhmin = hmin; 1364 r_8 dhmax = hmax; 1365 r_8 hb,hbmn,hbmx; 1352 1366 1353 1367 if(hdyn==0) hdyn = 100; … … 1355 1369 cout << "~Histo::Print " 1356 1370 << " nHist=" << nHist << " nEntries=" << nEntries 1357 << " under=" << under << " over=" << over << endl;1358 cout << " bins=" << bins1359 << " min=" << m in << " max=" << max1371 << " mUnder=" << mUnder << " mOver=" << mOver << endl; 1372 cout << " mBins=" << mBins 1373 << " min=" << mMin << " mMax=" << mMax 1360 1374 << " binWidth=" << binWidth << endl; 1361 1375 cout << " mean=" << Mean() << " r.m.s=" << Sigma() … … 1364 1378 if(hdyn<0 || pflag<0 ) return; 1365 1379 1366 if(dhmax<=dhmin) { if(hmin != 0.) dhmin = (double)VMin(); else dhmin=0.;1367 dhmax = (double)VMax(); }1380 if(dhmax<=dhmin) { if(hmin != 0.) dhmin = VMin(); else dhmin=0.; 1381 dhmax = VMax(); } 1368 1382 if(dhmin>dhmax) return; 1369 1383 if(dhmin==dhmax) {dhmin -= 1.; dhmax += 1.;} 1370 double wb = (dhmax-dhmin) / (double)hdyn;1384 r_8 wb = (dhmax-dhmin) / hdyn; 1371 1385 1372 1386 // determination de la position de la valeur zero 1373 int i0 = (int)(-dhmin/wb);1387 int_4 i0 = (int_4)(-dhmin/wb); 1374 1388 1375 1389 // si le zero est dans la dynamique, … … 1389 1403 dhmax -= hbmx; 1390 1404 } 1391 wb = (dhmax-dhmin) / (double)hdyn;1392 i0 = (int )(-dhmin/wb);1405 wb = (dhmax-dhmin) / hdyn; 1406 i0 = (int_4)(-dhmin/wb); 1393 1407 } 1394 1408 } … … 1401 1415 1402 1416 // premiere ligne 1403 {for(int i=0;i<hdyn;i++) s[i] = '=';}1417 {for(int_4 i=0;i<hdyn;i++) s[i] = '=';} 1404 1418 if( 0 <= i0 && i0 < hdyn ) s[i0] = '0'; 1405 1419 if(pflag&1) fprintf( fp, "===="); 1406 1420 fprintf( fp, "======================"); 1407 if(pflag&2 && err2!=NULL) fprintf( fp, "===========");1421 if(pflag&2 && mErr2!=NULL) fprintf( fp, "==========="); 1408 1422 fprintf( fp, "==%s\n",s); 1409 1423 1410 1424 // histogramme 1411 {for(int i=il;i<=ih;i++) {1412 for(int k=0;k<hdyn;k++) s[k] = ' ';1425 {for(int_4 i=il;i<=ih;i++) { 1426 for(int_4 k=0;k<hdyn;k++) s[k] = ' '; 1413 1427 hb = (*this)(i); 1414 1428 1415 1429 //choix du bin (hdyn bin entre [dhmin,dhmax[ 1416 int ii = (int)( (hb-dhmin)/wb );1430 int_4 ii = (int_4)( (hb-dhmin)/wb ); 1417 1431 if(ii<0) ii = 0; else if(ii>=hdyn) ii = hdyn-1; 1418 1432 … … 1424 1438 if(i0<0) { 1425 1439 // courbe entierement positive 1426 for(int k=0;k<=ii;k++) s[k] = 'X';1440 for(int_4 k=0;k<=ii;k++) s[k] = 'X'; 1427 1441 } else if(i0>=hdyn) { 1428 1442 // courbe entierement negative 1429 for(int k=hdyn-1;k>=ii;k--) s[k] = 'X';1443 for(int_4 k=hdyn-1;k>=ii;k--) s[k] = 'X'; 1430 1444 } else { 1431 1445 // courbe positive et negative 1432 1446 s[i0] = '|'; 1433 if(ii>i0) for(int k=i0+1;k<=ii;k++) s[k] = 'X';1434 else if(ii<i0) for(int k=ii;k<i0;k++) s[k] = 'X';1447 if(ii>i0) for(int_4 k=i0+1;k<=ii;k++) s[k] = 'X'; 1448 else if(ii<i0) for(int_4 k=ii;k<i0;k++) s[k] = 'X'; 1435 1449 else s[ii] = 'X'; 1436 1450 } 1437 1451 1438 1452 // valeur a mettre dans le bin le plus haut/bas 1439 int ib;1440 if(hb>0.) ib = (int )( 10.*(hb-hbmn)/(hbmx-hbmn) );1441 else if(hb<0.) ib = (int )( 10.*(hbmx-hb)/(hbmx-hbmn) );1453 int_4 ib; 1454 if(hb>0.) ib = (int_4)( 10.*(hb-hbmn)/(hbmx-hbmn) ); 1455 else if(hb<0.) ib = (int_4)( 10.*(hbmx-hb)/(hbmx-hbmn) ); 1442 1456 else ib = -1; 1443 1457 if(ib==-1) s[ii] = '|'; 1444 1458 else if(ib==0) s[ii] = '.'; 1445 else if(ib>0 && ib<10) s[ii] = (char)((int ) '0' + ib);1459 else if(ib>0 && ib<10) s[ii] = (char)((int_4) '0' + ib); 1446 1460 1447 1461 // traitement des saturations +/- … … 1451 1465 if(pflag&1) fprintf( fp, "%3d ",i); 1452 1466 fprintf( fp, "%10.4g %10.4g ",BinCenter(i),hb); 1453 if(pflag&2 && err2!=NULL) fprintf( fp, "%10.4g ",Error(i));1467 if(pflag&2 && mErr2!=NULL) fprintf( fp, "%10.4g ",Error(i)); 1454 1468 fprintf( fp, "= %s\n",s); 1455 1469 }} 1456 1470 1457 1471 // derniere ligne 1458 for(int i=0;i<hdyn;i++) s[i] = '=';1472 for(int_4 i=0;i<hdyn;i++) s[i] = '='; 1459 1473 if( 0 <= i0 && i0 < hdyn ) s[i0] = '0'; 1460 1474 if(pflag&1) fprintf( fp, "===="); 1461 1475 fprintf( fp, "======================"); 1462 if(pflag&2 && err2!=NULL) fprintf( fp, "===========");1476 if(pflag&2 && mErr2!=NULL) fprintf( fp, "==========="); 1463 1477 fprintf( fp, "==%s\n",s); 1464 1478 1465 1479 // valeur basse des bins (sur ["ndig-1" digits + signe] = ndig char (>=3)) 1466 const int ndig = 7;1480 const int_4 ndig = 7; 1467 1481 char sn[2*ndig+10]; 1468 1482 hb = ( fabs(dhmin) > fabs(dhmax) ) ? fabs(dhmin) : fabs(dhmax); 1469 int n;1470 if( hb > 0. ) n = (int )(log10(hb*1.00000001)); else n = 1;1471 double scaledig = pow(10.,(double) ndig-2);1472 double expo = scaledig/pow(10.,(double) n);1483 int_4 n; 1484 if( hb > 0. ) n = (int_4)(log10(hb*1.00000001)); else n = 1; 1485 r_8 scaledig = pow(10.,(r_8) ndig-2); 1486 r_8 expo = scaledig/pow(10.,(r_8) n); 1473 1487 // cout <<"n="<<n<<" hb="<<hb<<" scaledig="<<scaledig<<" expo="<<expo<<endl; 1474 for(int k=0;k<ndig;k++) {1475 for(int i=0;i<hdyn;i++) s[i] = ' ';1476 {for(int i=0;i<hdyn;i++) {1477 n = (int )( (dhmin + i*wb)*expo );1478 for(int j=0;j<2*ndig+10;j++) sn[j] = ' ';1488 for(int_4 k=0;k<ndig;k++) { 1489 for(int_4 i=0;i<hdyn;i++) s[i] = ' '; 1490 {for(int_4 i=0;i<hdyn;i++) { 1491 n = (int_4)( (dhmin + i*wb)*expo ); 1492 for(int_4 j=0;j<2*ndig+10;j++) sn[j] = ' '; 1479 1493 sprintf(sn,"%d%c",n,'\0'); 1480 1494 strip(sn,'B',' '); 1481 1495 // cout <<"n="<<n<<" sn=("<<sn<<") l="<<strlen(sn)<<" k="<<k; 1482 if( (int ) strlen(sn) > k ) s[i] = sn[k];1496 if( (int_4) strlen(sn) > k ) s[i] = sn[k]; 1483 1497 // cout <<" s=("<<s<<")"<<endl; 1484 1498 }} 1485 1499 if(pflag&1) fprintf( fp, " "); 1486 1500 fprintf( fp, " "); 1487 if(pflag&2 && err2!=NULL) fprintf( fp, " ");1501 if(pflag&2 && mErr2!=NULL) fprintf( fp, " "); 1488 1502 fprintf( fp, " %s\n",s); 1489 1503 } … … 1497 1511 Impression de l'histogramme sur stdout 1498 1512 */ 1499 void Histo::Print(int hdyn,float hmin, float hmax,intpflag,1500 int il, intih)1513 void Histo::Print(int_4 hdyn,r_8 hmin, r_8 hmax,int_4 pflag, 1514 int_4 il, int_4 ih) 1501 1515 { 1502 1516 Histo::PrintF(stdout, hdyn, hmin, hmax, pflag, il, ih); … … 1523 1537 1524 1538 // Lecture des valeurs 1525 is.Get(dobj-> bins);1539 is.Get(dobj->mBins); 1526 1540 is.Get(dobj->nEntries); 1527 1541 int_4 errok; … … 1529 1543 1530 1544 is.Get(dobj->binWidth); 1531 is.Get(dobj->m in);1532 is.Get(dobj->m ax);1533 is.Get(dobj-> under);1534 is.Get(dobj-> over);1545 is.Get(dobj->mMin); 1546 is.Get(dobj->mMax); 1547 is.Get(dobj->mUnder); 1548 is.Get(dobj->mOver); 1535 1549 1536 1550 is.Get(dobj->nHist); 1537 1551 1538 1552 // Lecture des donnees Histo 1D 1539 dobj-> data = new float[dobj->bins];1553 dobj->mData = new r_8[dobj->mBins]; 1540 1554 is.GetLine(strg, 255); 1541 is.Get(dobj-> data, dobj->bins);1555 is.Get(dobj->mData, dobj->mBins); 1542 1556 1543 1557 // Lecture des erreurs 1544 1558 if(errok) { 1545 1559 is.GetLine(strg, 255); 1546 dobj-> err2 = new double[dobj->bins];1547 is.Get(dobj-> err2, dobj->bins);1560 dobj->mErr2 = new r_8[dobj->mBins]; 1561 is.Get(dobj->mErr2, dobj->mBins); 1548 1562 } 1549 1563 … … 1557 1571 1558 1572 // Que faut-il ecrire? 1559 int_4 errok = (dobj-> err2) ? 1 : 0;1573 int_4 errok = (dobj->mErr2) ? 1 : 0; 1560 1574 1561 1575 // Ecriture entete pour identifier facilement 1562 sprintf(strg," bins=%6d NEnt=%15d errok=%1d",dobj->bins,dobj->nEntries,errok);1576 sprintf(strg,"mBins=%6d NEnt=%15d errok=%1d",dobj->mBins,dobj->nEntries,errok); 1563 1577 os.PutLine(strg); 1564 sprintf(strg,"binw=%g m in=%g max=%g",dobj->binWidth,dobj->min,dobj->max);1578 sprintf(strg,"binw=%g mMin=%g mMax=%g",dobj->binWidth,dobj->mMin,dobj->mMax); 1565 1579 os.PutLine(strg); 1566 sprintf(strg, " under=%g over=%g nHist=%g",dobj->under,dobj->over,dobj->nHist);1580 sprintf(strg, "mUnder=%g mOver=%g nHist=%g",dobj->mUnder,dobj->mOver,dobj->nHist); 1567 1581 os.PutLine(strg); 1568 1582 1569 1583 // Ecriture des valeurs 1570 os.Put(dobj-> bins);1584 os.Put(dobj->mBins); 1571 1585 os.Put(dobj->nEntries); 1572 1586 os.Put(errok); 1573 1587 1574 1588 os.Put(dobj->binWidth); 1575 os.Put(dobj->m in);1576 os.Put(dobj->m ax);1577 os.Put(dobj-> under);1578 os.Put(dobj-> over);1589 os.Put(dobj->mMin); 1590 os.Put(dobj->mMax); 1591 os.Put(dobj->mUnder); 1592 os.Put(dobj->mOver); 1579 1593 1580 1594 os.Put(dobj->nHist); 1581 1595 1582 1596 // Ecriture des donnees Histo 1D 1583 sprintf(strg,"Histo: Tableau des donnees %d",dobj-> bins);1597 sprintf(strg,"Histo: Tableau des donnees %d",dobj->mBins); 1584 1598 os.PutLine(strg); 1585 os.Put(dobj-> data, dobj->bins);1599 os.Put(dobj->mData, dobj->mBins); 1586 1600 1587 1601 // Ecriture des erreurs 1588 1602 if(errok) { 1589 sprintf(strg,"Histo: Tableau des erreurs %d",dobj-> bins);1603 sprintf(strg,"Histo: Tableau des erreurs %d",dobj->mBins); 1590 1604 os.PutLine(strg); 1591 os.Put(dobj-> err2, dobj->bins);1605 os.Put(dobj->mErr2, dobj->mBins); 1592 1606 } 1593 1607
Note:
See TracChangeset
for help on using the changeset viewer.