00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include <stdlib.h>
00013 #include <string.h>
00014 #include <stdio.h>
00015 #include <ctype.h>
00016
00017 #include "Riostream.h"
00018 #include "TROOT.h"
00019 #include "TClass.h"
00020 #include "TMath.h"
00021 #include "THashList.h"
00022 #include "TH1.h"
00023 #include "TH2.h"
00024 #include "TF2.h"
00025 #include "TF3.h"
00026 #include "TPluginManager.h"
00027 #include "TVirtualPad.h"
00028 #include "TRandom.h"
00029 #include "TVirtualFitter.h"
00030 #include "THLimitsFinder.h"
00031 #include "TProfile.h"
00032 #include "TStyle.h"
00033 #include "TVectorF.h"
00034 #include "TVectorD.h"
00035 #include "TBrowser.h"
00036 #include "TObjString.h"
00037 #include "TError.h"
00038 #include "TVirtualHistPainter.h"
00039 #include "TVirtualFFT.h"
00040 #include "TSystem.h"
00041
00042 #include "HFitInterface.h"
00043 #include "Fit/DataRange.h"
00044 #include "Math/MinimizerOptions.h"
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522 TF1 *gF1=0;
00523
00524 Int_t TH1::fgBufferSize = 1000;
00525 Bool_t TH1::fgAddDirectory = kTRUE;
00526 Bool_t TH1::fgDefaultSumw2 = kFALSE;
00527 Bool_t TH1::fgStatOverflows= kFALSE;
00528
00529 extern void H1InitGaus();
00530 extern void H1InitExpo();
00531 extern void H1InitPolynom();
00532 extern void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a);
00533 extern void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail);
00534 extern void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b);
00535
00536
00537 class DifferentNumberOfBins: public std::exception {};
00538 class DifferentAxisLimits: public std::exception {};
00539 class DifferentBinLimits: public std::exception {};
00540
00541 ClassImp(TH1)
00542
00543
00544 TH1::TH1(): TNamed(), TAttLine(), TAttFill(), TAttMarker()
00545 {
00546
00547
00548 fDirectory = 0;
00549 fFunctions = new TList;
00550 fNcells = 0;
00551 fIntegral = 0;
00552 fPainter = 0;
00553 fEntries = 0;
00554 fNormFactor = 0;
00555 fTsumw = fTsumw2=fTsumwx=fTsumwx2=0;
00556 fMaximum = -1111;
00557 fMinimum = -1111;
00558 fBufferSize = 0;
00559 fBuffer = 0;
00560 fXaxis.SetName("xaxis");
00561 fYaxis.SetName("yaxis");
00562 fZaxis.SetName("zaxis");
00563 fXaxis.SetParent(this);
00564 fYaxis.SetParent(this);
00565 fZaxis.SetParent(this);
00566 UseCurrentStyle();
00567 }
00568
00569
00570 TH1::~TH1()
00571 {
00572
00573
00574
00575 if (!TestBit(kNotDeleted)) {
00576 return;
00577 }
00578 delete[] fIntegral;
00579 fIntegral = 0;
00580 delete[] fBuffer;
00581 fBuffer = 0;
00582 if (fFunctions) {
00583 fFunctions->SetBit(kInvalidObject);
00584 TObject* obj = 0;
00585
00586
00587
00588
00589
00590
00591
00592 while ((obj = fFunctions->First())) {
00593 while(fFunctions->Remove(obj)) { }
00594 if (!obj->TestBit(kNotDeleted)) {
00595 break;
00596 }
00597 delete obj;
00598 obj = 0;
00599 }
00600 delete fFunctions;
00601 fFunctions = 0;
00602 }
00603 if (fDirectory) {
00604 fDirectory->Remove(this);
00605 fDirectory = 0;
00606 }
00607 delete fPainter;
00608 fPainter = 0;
00609 }
00610
00611
00612 TH1::TH1(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
00613 :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
00614 {
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 Build();
00637 if (nbins <= 0) nbins = 1;
00638 fXaxis.Set(nbins,xlow,xup);
00639 fNcells = fXaxis.GetNbins()+2;
00640 }
00641
00642
00643 TH1::TH1(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
00644 :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
00645 {
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660 Build();
00661 if (nbins <= 0) nbins = 1;
00662 if (xbins) fXaxis.Set(nbins,xbins);
00663 else fXaxis.Set(nbins,0,1);
00664 fNcells = fXaxis.GetNbins()+2;
00665 }
00666
00667
00668 TH1::TH1(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
00669 :TNamed(name,title), TAttLine(), TAttFill(), TAttMarker()
00670 {
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685 Build();
00686 if (nbins <= 0) nbins = 1;
00687 if (xbins) fXaxis.Set(nbins,xbins);
00688 else fXaxis.Set(nbins,0,1);
00689 fNcells = fXaxis.GetNbins()+2;
00690 }
00691
00692
00693 TH1::TH1(const TH1 &h) : TNamed(), TAttLine(), TAttFill(), TAttMarker()
00694 {
00695
00696
00697
00698 ((TH1&)h).Copy(*this);
00699 }
00700
00701
00702 Bool_t TH1::AddDirectoryStatus()
00703 {
00704
00705 return fgAddDirectory;
00706 }
00707
00708
00709 void TH1::Browse(TBrowser *b)
00710 {
00711
00712
00713 Draw(b ? b->GetDrawOption() : "");
00714 gPad->Update();
00715 }
00716
00717
00718
00719 void TH1::Build()
00720 {
00721
00722
00723
00724 fDirectory = 0;
00725 fPainter = 0;
00726 fIntegral = 0;
00727 fEntries = 0;
00728 fNormFactor = 0;
00729 fTsumw = fTsumw2=fTsumwx=fTsumwx2=0;
00730 fMaximum = -1111;
00731 fMinimum = -1111;
00732 fBufferSize = 0;
00733 fBuffer = 0;
00734 fXaxis.SetName("xaxis");
00735 fYaxis.SetName("yaxis");
00736 fZaxis.SetName("zaxis");
00737 fYaxis.Set(1,0.,1.);
00738 fZaxis.Set(1,0.,1.);
00739 fXaxis.SetParent(this);
00740 fYaxis.SetParent(this);
00741 fZaxis.SetParent(this);
00742
00743 SetTitle(fTitle.Data());
00744
00745 fFunctions = new TList;
00746
00747 UseCurrentStyle();
00748
00749 if (TH1::AddDirectoryStatus()) {
00750 fDirectory = gDirectory;
00751 if (fDirectory) {
00752 fDirectory->Append(this,kTRUE);
00753 }
00754 }
00755 }
00756
00757
00758 void TH1::Add(TF1 *f1, Double_t c1, Option_t *option)
00759 {
00760
00761
00762
00763
00764
00765
00766
00767
00768
00769
00770
00771
00772 if (!f1) {
00773 Error("Add","Attempt to add a non-existing function");
00774 return;
00775 }
00776
00777 TString opt = option;
00778 opt.ToLower();
00779 Bool_t integral = kFALSE;
00780 if (opt.Contains("i") && fDimension ==1) integral = kTRUE;
00781
00782 Int_t nbinsx = GetNbinsX();
00783 Int_t nbinsy = GetNbinsY();
00784 Int_t nbinsz = GetNbinsZ();
00785 if (fDimension < 2) nbinsy = -1;
00786 if (fDimension < 3) nbinsz = -1;
00787
00788
00789 Double_t s1[10];
00790 Int_t i;
00791 for (i=0;i<10;i++) {s1[i] = 0;}
00792 PutStats(s1);
00793 SetMinimum();
00794 SetMaximum();
00795
00796
00797 Int_t bin, binx, biny, binz;
00798 Double_t cu=0;
00799 Double_t xx[3];
00800 Double_t *params = 0;
00801 f1->InitArgs(xx,params);
00802 for (binz=0;binz<=nbinsz+1;binz++) {
00803 xx[2] = fZaxis.GetBinCenter(binz);
00804 for (biny=0;biny<=nbinsy+1;biny++) {
00805 xx[1] = fYaxis.GetBinCenter(biny);
00806 for (binx=0;binx<=nbinsx+1;binx++) {
00807 xx[0] = fXaxis.GetBinCenter(binx);
00808 if (!f1->IsInside(xx)) continue;
00809 TF1::RejectPoint(kFALSE);
00810 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
00811 if (integral) {
00812 xx[0] = fXaxis.GetBinLowEdge(binx);
00813 cu = c1*f1->EvalPar(xx);
00814 cu += c1*f1->Integral(fXaxis.GetBinLowEdge(binx),fXaxis.GetBinUpEdge(binx))*fXaxis.GetBinWidth(binx);
00815 } else {
00816 cu = c1*f1->EvalPar(xx);
00817 }
00818 if (TF1::RejectedPoint()) continue;
00819 Double_t error1 = GetBinError(bin);
00820 AddBinContent(bin,cu);
00821 if (fSumw2.fN) {
00822
00823 fSumw2.fArray[bin] = error1*error1;
00824 }
00825 }
00826 }
00827 }
00828 }
00829
00830
00831 void TH1::Add(const TH1 *h1, Double_t c1)
00832 {
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853 if (!h1) {
00854 Error("Add","Attempt to add a non-existing histogram");
00855 return;
00856 }
00857
00858 Int_t nbinsx = GetNbinsX();
00859 Int_t nbinsy = GetNbinsY();
00860 Int_t nbinsz = GetNbinsZ();
00861
00862 try {
00863 CheckConsistency(this,h1);
00864 } catch(DifferentNumberOfBins&) {
00865 Error("Add","Attempt to add histograms with different number of bins");
00866 return;
00867 } catch(DifferentAxisLimits&) {
00868 Warning("Add","Attempt to add histograms with different axis limits");
00869 } catch(DifferentBinLimits&) {
00870 Warning("Add","Attempt to add histograms with different bin limits");
00871 }
00872
00873 if (fDimension < 2) nbinsy = -1;
00874 if (fDimension < 3) nbinsz = -1;
00875
00876
00877 if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
00878
00879
00880 Double_t entries = TMath::Abs( GetEntries() + c1 * h1->GetEntries() );
00881 Double_t s1[kNstat], s2[kNstat];
00882
00883
00884 Bool_t resetStats = (c1 < 0);
00885 if (!resetStats) {
00886
00887
00888 for (Int_t i=0;i<kNstat;i++) {
00889 s1[i] = 0;
00890 s2[i] = 0;
00891 }
00892 GetStats(s1);
00893 h1->GetStats(s2);
00894 for (Int_t i=0;i<kNstat;i++) {
00895 if (i == 1) s1[i] += c1*c1*s2[i];
00896 else s1[i] += c1*s2[i];
00897 }
00898 }
00899
00900 SetMinimum();
00901 SetMaximum();
00902
00903
00904 Int_t bin, binx, biny, binz;
00905 Double_t cu;
00906 Double_t factor =1;
00907 if (h1->GetNormFactor() != 0) factor = h1->GetNormFactor()/h1->GetSumOfWeights();;
00908 for (binz=0;binz<=nbinsz+1;binz++) {
00909 for (biny=0;biny<=nbinsy+1;biny++) {
00910 for (binx=0;binx<=nbinsx+1;binx++) {
00911 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
00912
00913 if (this->TestBit(kIsAverage) && h1->TestBit(kIsAverage)) {
00914 Double_t y1 = h1->GetBinContent(bin);
00915 Double_t y2 = this->GetBinContent(bin);
00916 Double_t e1 = h1->GetBinError(bin);
00917 Double_t e2 = this->GetBinError(bin);
00918 Double_t w1 = 1., w2 = 1.;
00919 if (e1 > 0) w1 = 1./(e1*e1);
00920 if (e2 > 0) w2 = 1./(e2*e2);
00921 SetBinContent(bin, (w1*y1 + w2*y2)/(w1 + w2));
00922 if (fSumw2.fN) fSumw2.fArray[bin] = 1./(w1 + w2);
00923 } else {
00924 cu = c1*factor*h1->GetBinContent(bin);
00925 AddBinContent(bin,cu);
00926 if (fSumw2.fN) {
00927 Double_t e1 = factor*h1->GetBinError(bin);
00928 fSumw2.fArray[bin] += c1*c1*e1*e1;
00929 }
00930 }
00931 }
00932 }
00933 }
00934
00935
00936 if (resetStats) {
00937
00938 ResetStats();
00939 }
00940 else {
00941 PutStats(s1);
00942 SetEntries(entries);
00943 }
00944 }
00945
00946
00947 void TH1::Add(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2)
00948 {
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968 if (!h1 || !h2) {
00969 Error("Add","Attempt to add a non-existing histogram");
00970 return;
00971 }
00972
00973 Bool_t normWidth = kFALSE;
00974 if (h1 == h2 && c2 < 0) {c2 = 0; normWidth = kTRUE;}
00975 Int_t nbinsx = GetNbinsX();
00976 Int_t nbinsy = GetNbinsY();
00977 Int_t nbinsz = GetNbinsZ();
00978
00979 try {
00980 CheckConsistency(h1,h2);
00981 CheckConsistency(this,h1);
00982 } catch(DifferentNumberOfBins&) {
00983 Error("Add","Attempt to add histograms with different number of bins");
00984 return;
00985 } catch(DifferentAxisLimits&) {
00986 Warning("Add","Attempt to add histograms with different axis limits");
00987 } catch(DifferentBinLimits&) {
00988 Warning("Add","Attempt to add histograms with different bin limits");
00989 }
00990
00991 if (fDimension < 2) nbinsy = -1;
00992 if (fDimension < 3) nbinsz = -1;
00993 if (fDimension < 3) nbinsz = -1;
00994
00995
00996 if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
00997
00998
00999 Double_t nEntries = TMath::Abs( c1*h1->GetEntries() + c2*h2->GetEntries() );
01000 Double_t s1[kNstat], s2[kNstat], s3[kNstat];
01001
01002
01003
01004 Bool_t resetStats = (c1*c2 < 0);
01005 if (!resetStats) {
01006
01007
01008 for (Int_t i=0;i<kNstat;i++) {
01009 s1[i] = 0;
01010 s2[i] = 0;
01011 }
01012 h1->GetStats(s1);
01013 h2->GetStats(s2);
01014 for (Int_t i=0;i<kNstat;i++) {
01015 if (i == 1) s3[i] = c1*c1*s1[i] + c2*c2*s2[i];
01016 else s3[i] = TMath::Abs(c1)*s1[i] + TMath::Abs(c2)*s2[i];
01017 }
01018 }
01019
01020 SetMinimum();
01021 SetMaximum();
01022
01023
01024
01025 ResetBit(kCanRebin);
01026
01027
01028
01029 Int_t bin, binx, biny, binz;
01030 Double_t cu;
01031 for (binz=0;binz<=nbinsz+1;binz++) {
01032 Double_t wz = h1->GetZaxis()->GetBinWidth(binz);
01033 for (biny=0;biny<=nbinsy+1;biny++) {
01034 Double_t wy = h1->GetYaxis()->GetBinWidth(biny);
01035 for (binx=0;binx<=nbinsx+1;binx++) {
01036 Double_t wx = h1->GetXaxis()->GetBinWidth(binx);
01037 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
01038
01039 if (h1->TestBit(kIsAverage) && h2->TestBit(kIsAverage)) {
01040 Double_t y1 = h1->GetBinContent(bin);
01041 Double_t y2 = h2->GetBinContent(bin);
01042 Double_t e1 = h1->GetBinError(bin);
01043 Double_t e2 = h2->GetBinError(bin);
01044 Double_t w1 = 1., w2 = 1.;
01045 if (e1 > 0) w1 = 1./(e1*e1);
01046 if (e2 > 0) w2 = 1./(e2*e2);
01047 SetBinContent(bin, (w1*y1 + w2*y2)/(w1 + w2));
01048 if (fSumw2.fN) fSumw2.fArray[bin] = 1./(w1 + w2);
01049 } else {
01050 if (normWidth) {
01051 Double_t w = wx*wy*wz;
01052 cu = c1*h1->GetBinContent(bin)/w;
01053 SetBinContent(bin,cu);
01054 if (fSumw2.fN) {
01055 Double_t e1 = h1->GetBinError(bin)/w;
01056 fSumw2.fArray[bin] = c1*c1*e1*e1;
01057 }
01058 } else {
01059 cu = c1*h1->GetBinContent(bin)+ c2*h2->GetBinContent(bin);
01060 SetBinContent(bin,cu);
01061 if (fSumw2.fN) {
01062 Double_t e1 = h1->GetBinError(bin);
01063 Double_t e2 = h2->GetBinError(bin);
01064 fSumw2.fArray[bin] = c1*c1*e1*e1 + c2*c2*e2*e2;
01065 }
01066 }
01067 }
01068 }
01069 }
01070 }
01071 if (resetStats) {
01072
01073 ResetStats();
01074 }
01075 else {
01076
01077 PutStats(s3);
01078 SetEntries(nEntries);
01079 }
01080 }
01081
01082
01083
01084 void TH1::AddBinContent(Int_t)
01085 {
01086
01087
01088 AbstractMethod("AddBinContent");
01089 }
01090
01091
01092 void TH1::AddBinContent(Int_t, Double_t)
01093 {
01094
01095
01096 AbstractMethod("AddBinContent");
01097 }
01098
01099
01100 void TH1::AddDirectory(Bool_t add)
01101 {
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113 fgAddDirectory = add;
01114 }
01115
01116
01117
01118 Int_t TH1::BufferEmpty(Int_t action)
01119 {
01120
01121
01122
01123
01124
01125
01126
01127
01128 if (!fBuffer) return 0;
01129 Int_t nbentries = (Int_t)fBuffer[0];
01130 if (!nbentries) return 0;
01131 Double_t *buffer = fBuffer;
01132 if (nbentries < 0) {
01133 if (action == 0) return 0;
01134 nbentries = -nbentries;
01135 fBuffer=0;
01136 Reset("ICE");
01137 fBuffer = buffer;
01138 }
01139 if (TestBit(kCanRebin) || (fXaxis.GetXmax() <= fXaxis.GetXmin())) {
01140
01141 Double_t xmin = fBuffer[2];
01142 Double_t xmax = xmin;
01143 for (Int_t i=1;i<nbentries;i++) {
01144 Double_t x = fBuffer[2*i+2];
01145 if (x < xmin) xmin = x;
01146 if (x > xmax) xmax = x;
01147 }
01148 if (fXaxis.GetXmax() <= fXaxis.GetXmin()) {
01149 THLimitsFinder::GetLimitsFinder()->FindGoodLimits(this,xmin,xmax);
01150 } else {
01151 fBuffer = 0;
01152 Int_t keep = fBufferSize; fBufferSize = 0;
01153 if (xmin < fXaxis.GetXmin()) RebinAxis(xmin,&fXaxis);
01154 if (xmax >= fXaxis.GetXmax()) RebinAxis(xmax,&fXaxis);
01155 fBuffer = buffer;
01156 fBufferSize = keep;
01157 }
01158 }
01159
01160 FillN(nbentries,&fBuffer[2],&fBuffer[1],2);
01161
01162 if (action > 0) { delete [] fBuffer; fBuffer = 0; fBufferSize = 0;}
01163 else {
01164 if (nbentries == (Int_t)fEntries) fBuffer[0] = -nbentries;
01165 else fBuffer[0] = 0;
01166 }
01167 return nbentries;
01168 }
01169
01170
01171 Int_t TH1::BufferFill(Double_t x, Double_t w)
01172 {
01173
01174
01175
01176
01177
01178 if (!fBuffer) return -2;
01179 Int_t nbentries = (Int_t)fBuffer[0];
01180 if (nbentries < 0) {
01181 nbentries = -nbentries;
01182 fBuffer[0] = nbentries;
01183 if (fEntries > 0) {
01184 Double_t *buffer = fBuffer; fBuffer=0;
01185 Reset();
01186 fBuffer = buffer;
01187 }
01188 }
01189 if (2*nbentries+2 >= fBufferSize) {
01190 BufferEmpty(1);
01191 return Fill(x,w);
01192 }
01193 fBuffer[2*nbentries+1] = w;
01194 fBuffer[2*nbentries+2] = x;
01195 fBuffer[0] += 1;
01196 return -2;
01197 }
01198
01199 bool CheckBinLimits(const TArrayD* h1Array, const TArrayD* h2Array)
01200 {
01201 Int_t fN = h1Array->fN;
01202 if ( fN != 0 ) {
01203 if ( h2Array->fN != fN ) {
01204 throw DifferentBinLimits();
01205 return false;
01206 }
01207 else {
01208 for ( int i = 0; i < fN; ++i ) {
01209 if ( ! TMath::AreEqualAbs( h1Array->GetAt(i), h2Array->GetAt(i), 1E-10 ) ) {
01210 throw DifferentBinLimits();
01211 return false;
01212 }
01213 }
01214 }
01215 }
01216
01217 return true;
01218 }
01219
01220
01221 bool TH1::CheckConsistency(const TH1* h1, const TH1* h2)
01222 {
01223
01224
01225 Int_t nbinsx = h1->GetNbinsX();
01226 Int_t nbinsy = h1->GetNbinsY();
01227 Int_t nbinsz = h1->GetNbinsZ();
01228
01229
01230 if (nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
01231 throw DifferentNumberOfBins();
01232 return false;
01233 }
01234
01235 if (h1->fXaxis.GetXmin() != h2->fXaxis.GetXmin() ||
01236 h1->fXaxis.GetXmax() != h2->fXaxis.GetXmax() ||
01237 h1->fYaxis.GetXmin() != h2->fYaxis.GetXmin() ||
01238 h1->fYaxis.GetXmax() != h2->fYaxis.GetXmax() ||
01239 h1->fZaxis.GetXmin() != h2->fZaxis.GetXmin() ||
01240 h1->fZaxis.GetXmax() != h2->fZaxis.GetXmax()) {
01241 throw DifferentAxisLimits();
01242 return false;
01243 }
01244
01245 bool ret = true;
01246
01247 ret &= CheckBinLimits(h1->GetXaxis()->GetXbins(), h2->GetXaxis()->GetXbins());
01248 ret &= CheckBinLimits(h1->GetYaxis()->GetXbins(), h2->GetYaxis()->GetXbins());
01249 ret &= CheckBinLimits(h1->GetZaxis()->GetXbins(), h2->GetZaxis()->GetXbins());
01250
01251 return ret;
01252 }
01253
01254
01255 Double_t TH1::Chi2Test(const TH1* h2, Option_t *option, Double_t *res) const
01256 {
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448
01449
01450
01451
01452
01453
01454
01455
01456
01457
01458
01459
01460
01461
01462
01463
01464
01465
01466
01467
01468
01469
01470
01471
01472
01473
01474
01475
01476
01477
01478
01479
01480
01481
01482
01483
01484
01485
01486
01487
01488
01489
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545 Double_t chi2 = 0;
01546 Int_t ndf = 0, igood = 0;
01547
01548 TString opt = option;
01549 opt.ToUpper();
01550
01551 Double_t prob = Chi2TestX(h2,chi2,ndf,igood,option,res);
01552
01553 if(opt.Contains("P")) {
01554 printf("Chi2 = %f, Prob = %g, NDF = %d, igood = %d\n", chi2,prob,ndf,igood);
01555 }
01556 if(opt.Contains("CHI2/NDF")) {
01557 if (ndf == 0) return 0;
01558 return chi2/ndf;
01559 }
01560 if(opt.Contains("CHI2")) {
01561 return chi2;
01562 }
01563
01564 return prob;
01565 }
01566
01567
01568 Double_t TH1::Chi2TestX(const TH1* h2, Double_t &chi2, Int_t &ndf, Int_t &igood, Option_t *option, Double_t *res) const
01569 {
01570
01571
01572
01573
01574
01575
01576
01577
01578
01579
01580
01581
01582
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592
01593
01594
01595
01596
01597
01598
01599
01600
01601
01602
01603
01604
01605
01606
01607
01608
01609
01610 Int_t i = 0, j=0, k = 0;
01611 Int_t i_start, i_end;
01612 Int_t j_start, j_end;
01613 Int_t k_start, k_end;
01614
01615 Double_t bin1, bin2;
01616 Double_t err1,err2;
01617 Double_t sum1=0, sum2=0;
01618 Double_t sumw1=0, sumw2=0;
01619
01620
01621 chi2 = 0;
01622 ndf = 0;
01623
01624 TString opt = option;
01625 opt.ToUpper();
01626
01627 TAxis *xaxis1 = this->GetXaxis();
01628 TAxis *xaxis2 = h2->GetXaxis();
01629 TAxis *yaxis1 = this->GetYaxis();
01630 TAxis *yaxis2 = h2->GetYaxis();
01631 TAxis *zaxis1 = this->GetZaxis();
01632 TAxis *zaxis2 = h2->GetZaxis();
01633
01634 Int_t nbinx1 = xaxis1->GetNbins();
01635 Int_t nbinx2 = xaxis2->GetNbins();
01636 Int_t nbiny1 = yaxis1->GetNbins();
01637 Int_t nbiny2 = yaxis2->GetNbins();
01638 Int_t nbinz1 = zaxis1->GetNbins();
01639 Int_t nbinz2 = zaxis2->GetNbins();
01640
01641
01642 if (this->GetDimension() != h2->GetDimension() ){
01643 Error("Chi2TestX","Histograms have different dimensions.");
01644 return 0;
01645 }
01646
01647
01648 if (nbinx1 != nbinx2) {
01649 Error("Chi2TestX","different number of x channels");
01650 }
01651 if (nbiny1 != nbiny2) {
01652 Error("Chi2TestX","different number of y channels");
01653 }
01654 if (nbinz1 != nbinz2) {
01655 Error("Chi2TestX","different number of z channels");
01656 }
01657
01658
01659 i_start = j_start = k_start = 1;
01660 i_end = nbinx1;
01661 j_end = nbiny1;
01662 k_end = nbinz1;
01663
01664 if (xaxis1->TestBit(TAxis::kAxisRange)) {
01665 i_start = xaxis1->GetFirst();
01666 i_end = xaxis1->GetLast();
01667 }
01668 if (yaxis1->TestBit(TAxis::kAxisRange)) {
01669 j_start = yaxis1->GetFirst();
01670 j_end = yaxis1->GetLast();
01671 }
01672 if (zaxis1->TestBit(TAxis::kAxisRange)) {
01673 k_start = zaxis1->GetFirst();
01674 k_end = zaxis1->GetLast();
01675 }
01676
01677
01678 if (opt.Contains("OF")) {
01679 if (this->GetDimension() == 3) k_end = ++nbinz1;
01680 if (this->GetDimension() >= 2) j_end = ++nbiny1;
01681 if (this->GetDimension() >= 1) i_end = ++nbinx1;
01682 }
01683
01684 if (opt.Contains("UF")) {
01685 if (this->GetDimension() == 3) k_start = 0;
01686 if (this->GetDimension() >= 2) j_start = 0;
01687 if (this->GetDimension() >= 1) i_start = 0;
01688 }
01689
01690 ndf = (i_end - i_start + 1)*(j_end - j_start + 1)*(k_end - k_start + 1) - 1;
01691
01692 Bool_t comparisonUU = opt.Contains("UU");
01693 Bool_t comparisonUW = opt.Contains("UW");
01694 Bool_t comparisonWW = opt.Contains("WW");
01695 Bool_t scaledHistogram = opt.Contains("NORM");
01696 if (scaledHistogram && !comparisonUU) {
01697 Info("Chi2TestX","NORM option should be used together with UU option. It is ignored");
01698 }
01699
01700 Stat_t s[kNstat];
01701 this->GetStats(s);
01702 double sumBinContent1 = s[0];
01703 double effEntries1 = (s[1] ? s[0]*s[0]/s[1] : 0.);
01704
01705 h2->GetStats(s);
01706 double sumBinContent2 = s[0];
01707 double effEntries2 = (s[1] ? s[0]*s[0]/s[1] : 0.);
01708
01709 if (!comparisonUU && !comparisonUW && !comparisonWW ) {
01710
01711 if (TMath::Abs(sumBinContent1 - effEntries1) < 1) {
01712 if ( TMath::Abs(sumBinContent2 - effEntries2) < 1) comparisonUU = true;
01713 else comparisonUW = true;
01714 }
01715 else comparisonWW = true;
01716 }
01717
01718 if (comparisonUW) {
01719 if (TMath::Abs(sumBinContent1 - effEntries1) >= 1) {
01720 Warning("Chi2TestX","First histogram is not unweighted and option UW has been requested");
01721 }
01722 }
01723 if ( (!scaledHistogram && comparisonUU) ) {
01724 if ( ( TMath::Abs(sumBinContent1 - effEntries1) >= 1) || (TMath::Abs(sumBinContent2 - effEntries2) >= 1) ) {
01725 Warning("Chi2TestX","Both histograms are not unweighted and option UU has been requested");
01726 }
01727 }
01728
01729
01730
01731 if (comparisonUU && scaledHistogram) {
01732 for (i=i_start; i<=i_end; i++) {
01733 for (j=j_start; j<=j_end; j++) {
01734 for (k=k_start; k<=k_end; k++) {
01735 bin1 = this->GetBinContent(i,j,k);
01736 bin2 = h2->GetBinContent(i,j,k);
01737 err1 = this->GetBinError(i,j,k);
01738 err2 = h2->GetBinError(i,j,k);
01739 if (err1 > 0 ) {
01740 bin1 *= bin1/(err1*err1);
01741
01742 bin1 = TMath::Floor(bin1+0.5);
01743 }
01744 else
01745 bin1 = 0;
01746
01747 if (err2 > 0) {
01748 bin2 *= bin2/(err2*err2);
01749
01750 bin2 = TMath::Floor(bin2+0.5);
01751 }
01752 else
01753 bin2 = 0;
01754
01755
01756 sum1 += bin1;
01757 sum2 += bin2;
01758 sumw1 += err1*err1;
01759 sumw2 += err2*err2;
01760 }
01761 }
01762 }
01763 if (sumw1 <= 0 || sumw2 <= 0) {
01764 Error("Chi2TestX","Cannot use option NORM when one histogram has all zero errors");
01765 return 0;
01766 }
01767
01768 } else {
01769 for (i=i_start; i<=i_end; i++) {
01770 for (j=j_start; j<=j_end; j++) {
01771 for (k=k_start; k<=k_end; k++) {
01772 sum1 += this->GetBinContent(i,j,k);
01773 sum2 += h2->GetBinContent(i,j,k);
01774 if ( comparisonWW ) {
01775 err1 = this->GetBinError(i,j,k);
01776 sumw1 += err1*err1;
01777 }
01778 if ( comparisonUW || comparisonWW ) {
01779 err2 = h2->GetBinError(i,j,k);
01780 sumw2 += err2*err2;
01781 }
01782 }
01783 }
01784 }
01785 }
01786
01787 if (sum1 == 0 || sum2 == 0) {
01788 Error("Chi2TestX","one histogram is empty");
01789 return 0;
01790 }
01791
01792 if ( comparisonWW && ( sumw1 <= 0 && sumw2 <=0 ) ){
01793 Error("Chi2TestX","Hist1 and Hist2 have both all zero errors\n");
01794 return 0;
01795 }
01796
01797
01798 Int_t m=0, n=0;
01799
01800
01801 if (comparisonUU) {
01802 Double_t sum = sum1 + sum2;
01803 for (i=i_start; i<=i_end; i++) {
01804 for (j=j_start; j<=j_end; j++) {
01805 for (k=k_start; k<=k_end; k++) {
01806 bin1 = this->GetBinContent(i,j,k);
01807 bin2 = h2->GetBinContent(i,j,k);
01808
01809 if ( (bin1 == 0) && (bin2 == 0) ) {
01810 --ndf;
01811 } else {
01812 if (scaledHistogram) {
01813
01814 err1 = this->GetBinError(i,j,k);
01815 if (err1 > 0 ) {
01816 bin1 *= bin1/(err1*err1);
01817
01818 bin1 = TMath::Floor(bin1+0.5);
01819 }
01820 else
01821 bin1 = 0;
01822
01823 err2 = h2->GetBinError(i,j,k);
01824 if (err2 > 0) {
01825 bin2 *= bin2/(err2*err2);
01826
01827 bin2 = TMath::Floor(bin2+0.5);
01828 }
01829 else
01830 bin2 = 0;
01831
01832 }
01833
01834
01835
01836 Double_t binsum = bin1 + bin2;
01837 Double_t nexp1 = binsum*sum1/sum;
01838
01839
01840
01841
01842 if (res)
01843 res[i-i_start] = (bin1-nexp1)/TMath::Sqrt(nexp1);
01844
01845 if (bin1 < 1) {
01846 m++;
01847 }
01848 if (bin2 < 1) {
01849 n++;
01850 }
01851
01852
01853 Double_t correc = (1-sum1/sum)*(1-binsum/sum);
01854 if (res) {
01855 res[i-i_start] /= TMath::Sqrt(correc);
01856 }
01857
01858 Double_t delta = sum2*bin1-sum1*bin2;
01859 chi2 += delta*delta/binsum;
01860
01861 }
01862 }
01863 }
01864 }
01865
01866 chi2 /= (sum1*sum2);
01867
01868 if (m) {
01869 igood += 1;
01870 Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
01871 }
01872 if (n) {
01873 igood += 2;
01874 Info("Chi2TestX","There is a bin in h2 with less than 1 event.\n");
01875 }
01876
01877 Double_t prob = TMath::Prob(chi2,ndf);
01878 return prob;
01879
01880 }
01881
01882
01883
01884 if ( comparisonUW ) {
01885 for (i=i_start; i<=i_end; i++) {
01886 for (j=j_start; j<=j_end; j++) {
01887 for (k=k_start; k<=k_end; k++) {
01888 Int_t x=0, y=0;
01889 bin1 = this->GetBinContent(i,j,k);
01890 bin2 = h2->GetBinContent(i,j,k);
01891 err2 = h2->GetBinError(i,j,k);
01892
01893 err2 *= err2;
01894
01895
01896 if ( (bin1 == 0) && (bin2 == 0) ) {
01897 --ndf;
01898 continue;
01899 }
01900
01901
01902 if (bin2 == 0 && err2 == 0) {
01903 if (sumw2 > 0) {
01904
01905
01906 err2 = sumw2/sum2;
01907 }
01908 else {
01909
01910
01911 Error("Chi2TestX","Hist2 has in bin %d,%d,%d zero content and all zero errors\n", i,j,k);
01912 chi2 = 0; return 0;
01913 }
01914 }
01915
01916
01917
01918 if (bin1 < 1) m++;
01919 if (err2 > 0 && bin2*bin2/err2 < 10) n++;
01920
01921 Double_t var1 = sum2*bin2 - sum1*err2;
01922 Double_t var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
01923
01924
01925 while (var1*var1+bin1 == 0 || var1+var2 == 0) {
01926 sum1++;
01927 bin1++;
01928 x++;
01929 y=1;
01930 var1 = sum2*bin2 - sum1*err2;
01931 var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
01932 }
01933 var2 = TMath::Sqrt(var2);
01934 while (var1+var2 == 0) {
01935 sum1++;
01936 bin1++;
01937 x++;
01938 y=1;
01939 var1 = sum2*bin2 - sum1*err2;
01940 var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
01941 while (var1*var1+bin1 == 0 || var1+var2 == 0) {
01942 sum1++;
01943 bin1++;
01944 x++;
01945 y=1;
01946 var1 = sum2*bin2 - sum1*err2;
01947 var2 = var1*var1 + 4*sum2*sum2*bin1*err2;
01948 }
01949 var2 = TMath::Sqrt(var2);
01950 }
01951
01952 Double_t probb = (var1+var2)/(2*sum2*sum2);
01953
01954 Double_t nexp1 = probb * sum1;
01955 Double_t nexp2 = probb * sum2;
01956
01957
01958
01959 Double_t delta1 = bin1 - nexp1;
01960 Double_t delta2 = bin2 - nexp2;
01961
01962 chi2 += delta1*delta1/nexp1;
01963
01964 if (err2 > 0) {
01965 chi2 += delta2*delta2/err2;
01966 }
01967
01968 if (res) {
01969 if (err2 > 0) {
01970 Double_t temp1 = sum2*err2/var2;
01971 Double_t temp2 = 1 + (sum1*err2 - sum2*bin2)/var2;
01972 temp2 = temp1*temp1*sum1*probb*(1-probb) + temp2*temp2*err2/4;
01973
01974 res[i-i_start] = - delta2/TMath::Sqrt(temp2);
01975 }
01976 else
01977 res[i-i_start] = delta1/TMath::Sqrt(nexp1);
01978
01979 }
01980 }
01981 }
01982 }
01983
01984 if (m) {
01985 igood += 1;
01986 Info("Chi2TestX","There is a bin in h1 with less than 1 event.\n");
01987 }
01988 if (n) {
01989 igood += 2;
01990 Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
01991 }
01992
01993 Double_t prob = TMath::Prob(chi2,ndf);
01994
01995 return prob;
01996 }
01997
01998
01999 if (comparisonWW) {
02000 for (i=i_start; i<=i_end; i++) {
02001 for (j=j_start; j<=j_end; j++) {
02002 for (k=k_start; k<=k_end; k++) {
02003 bin1 = this->GetBinContent(i,j,k);
02004 bin2 = h2->GetBinContent(i,j,k);
02005 err1 = this->GetBinError(i,j,k);
02006 err2 = h2->GetBinError(i,j,k);
02007 err1 *= err1;
02008 err2 *= err2;
02009
02010
02011 if ( (bin1 == 0) && (bin2 == 0) ) {
02012 --ndf;
02013 continue;
02014 }
02015 if ( (err1 == 0) && (err2 == 0) ) {
02016
02017 Error("Chi2TestX","h1 and h2 both have bin %d,%d,%d with all zero errors\n", i,j,k);
02018 chi2 = 0; return 0;
02019 }
02020
02021 Double_t sigma = sum1*sum1*err2 + sum2*sum2*err1;
02022 Double_t delta = sum2*bin1 - sum1*bin2;
02023 chi2 += delta*delta/sigma;
02024
02025
02026
02027 if (res) {
02028 Double_t temp = bin1*sum1*err2 + bin2*sum2*err1;
02029 Double_t probb = temp/sigma;
02030 Double_t z = 0;
02031 if (err1 > err2 ) {
02032 Double_t d1 = (bin1 - sum1 * probb);
02033 Double_t s1 = err1* ( 1. - err2 * sum1 * sum1 / sigma );
02034 z = d1/ TMath::Sqrt(s1);
02035 }
02036 else {
02037 Double_t d2 = (bin2 - sum2 * probb);
02038 Double_t s2 = err2* ( 1. - err1 * sum2 * sum2 / sigma );
02039 z = -d2/ TMath::Sqrt(s2);
02040 }
02041
02042 res[i-i_start] = z;
02043 }
02044
02045 if (err1 > 0 && bin1*bin1/err1 < 10) m++;
02046 if (err2 > 0 && bin2*bin2/err2 < 10) n++;
02047 }
02048 }
02049 }
02050 if (m) {
02051 igood += 1;
02052 Info("Chi2TestX","There is a bin in h1 with less than 10 effective events.\n");
02053 }
02054 if (n) {
02055 igood += 2;
02056 Info("Chi2TestX","There is a bin in h2 with less than 10 effective events.\n");
02057 }
02058 Double_t prob = TMath::Prob(chi2,ndf);
02059 return prob;
02060 }
02061 return 0;
02062 }
02063
02064
02065 Double_t TH1::ComputeIntegral()
02066 {
02067
02068
02069
02070
02071
02072
02073
02074 Int_t bin, binx, biny, binz, ibin;
02075
02076
02077 if (fIntegral) delete [] fIntegral;
02078
02079
02080 Int_t nbinsx = GetNbinsX();
02081 Int_t nbinsy = GetNbinsY();
02082 Int_t nbinsz = GetNbinsZ();
02083 Int_t nxy = nbinsx*nbinsy;
02084 Int_t nbins = nxy*nbinsz;
02085
02086 fIntegral = new Double_t[nbins+2];
02087 ibin = 0;
02088 fIntegral[ibin] = 0;
02089 for (binz=1;binz<=nbinsz;binz++) {
02090 for (biny=1;biny<=nbinsy;biny++) {
02091 for (binx=1;binx<=nbinsx;binx++) {
02092 ibin++;
02093 bin = GetBin(binx, biny, binz);
02094 fIntegral[ibin] = fIntegral[ibin-1] + GetBinContent(bin);
02095 }
02096 }
02097 }
02098
02099
02100 if (fIntegral[nbins] == 0 ) {
02101 Error("ComputeIntegral", "Integral = zero"); return 0;
02102 }
02103 for (bin=1;bin<=nbins;bin++) fIntegral[bin] /= fIntegral[nbins];
02104 fIntegral[nbins+1] = fEntries;
02105 return fIntegral[nbins];
02106 }
02107
02108
02109 Double_t *TH1::GetIntegral()
02110 {
02111
02112
02113
02114 if (!fIntegral) ComputeIntegral();
02115 return fIntegral;
02116 }
02117
02118
02119 void TH1::Copy(TObject &obj) const
02120 {
02121
02122
02123
02124
02125
02126
02127 if (((TH1&)obj).fDirectory) {
02128
02129
02130
02131 ((TH1&)obj).fDirectory->Remove(&obj);
02132 ((TH1&)obj).fDirectory = 0;
02133 }
02134 TNamed::Copy(obj);
02135 ((TH1&)obj).fDimension = fDimension;
02136 ((TH1&)obj).fNormFactor= fNormFactor;
02137 ((TH1&)obj).fNcells = fNcells;
02138 ((TH1&)obj).fBarOffset = fBarOffset;
02139 ((TH1&)obj).fBarWidth = fBarWidth;
02140 ((TH1&)obj).fOption = fOption;
02141 ((TH1&)obj).fBuffer = 0;
02142 ((TH1&)obj).fBufferSize= fBufferSize;
02143 Int_t i;
02144 if (fBuffer) {
02145 Double_t *buf = new Double_t[fBufferSize];
02146 for (i=0;i<fBufferSize;i++) buf[i] = fBuffer[i];
02147 ((TH1&)obj).fBuffer = buf;
02148 }
02149
02150 TArray* a = dynamic_cast<TArray*>(&obj);
02151 if (a) a->Set(fNcells);
02152 Int_t canRebin = ((TH1&)obj).TestBit(kCanRebin);
02153 ((TH1&)obj).ResetBit(kCanRebin);
02154 for (i=0;i<fNcells;i++) ((TH1&)obj).SetBinContent(i,this->GetBinContent(i));
02155 if (canRebin) ((TH1&)obj).SetBit(kCanRebin);
02156 ((TH1&)obj).fEntries = fEntries;
02157
02158 ((TH1&)obj).fTsumw = fTsumw;
02159 ((TH1&)obj).fTsumw2 = fTsumw2;
02160 ((TH1&)obj).fTsumwx = fTsumwx;
02161 ((TH1&)obj).fTsumwx2 = fTsumwx2;
02162 ((TH1&)obj).fMaximum = fMaximum;
02163 ((TH1&)obj).fMinimum = fMinimum;
02164
02165 TAttLine::Copy(((TH1&)obj));
02166 TAttFill::Copy(((TH1&)obj));
02167 TAttMarker::Copy(((TH1&)obj));
02168 fXaxis.Copy(((TH1&)obj).fXaxis);
02169 fYaxis.Copy(((TH1&)obj).fYaxis);
02170 fZaxis.Copy(((TH1&)obj).fZaxis);
02171 ((TH1&)obj).fXaxis.SetParent(&obj);
02172 ((TH1&)obj).fYaxis.SetParent(&obj);
02173 ((TH1&)obj).fZaxis.SetParent(&obj);
02174 fContour.Copy(((TH1&)obj).fContour);
02175 fSumw2.Copy(((TH1&)obj).fSumw2);
02176
02177 if (fgAddDirectory && gDirectory) {
02178 gDirectory->Append(&obj);
02179 ((TH1&)obj).fDirectory = gDirectory;
02180 }
02181 }
02182
02183
02184 void TH1::DirectoryAutoAdd(TDirectory *dir)
02185 {
02186
02187
02188
02189
02190
02191
02192
02193 Bool_t addStatus = TH1::AddDirectoryStatus();
02194 if (addStatus) {
02195 SetDirectory(dir);
02196 if (dir) {
02197 ResetBit(kCanDelete);
02198 }
02199 }
02200 }
02201
02202
02203 Int_t TH1::DistancetoPrimitive(Int_t px, Int_t py)
02204 {
02205
02206
02207
02208
02209
02210
02211
02212
02213
02214
02215
02216
02217 if (!fPainter) return 9999;
02218 return fPainter->DistancetoPrimitive(px,py);
02219 }
02220
02221
02222 void TH1::Divide(TF1 *f1, Double_t c1)
02223 {
02224
02225
02226
02227
02228
02229
02230
02231
02232 if (!f1) {
02233 Error("Add","Attempt to divide by a non-existing function");
02234 return;
02235 }
02236
02237 Int_t nbinsx = GetNbinsX();
02238 Int_t nbinsy = GetNbinsY();
02239 Int_t nbinsz = GetNbinsZ();
02240 if (fDimension < 2) nbinsy = -1;
02241 if (fDimension < 3) nbinsz = -1;
02242
02243
02244 SetMinimum();
02245 SetMaximum();
02246
02247
02248
02249 ResetBit(kCanRebin);
02250
02251
02252
02253 Int_t bin, binx, biny, binz;
02254 Double_t cu,w;
02255 Double_t xx[3];
02256 Double_t *params = 0;
02257 f1->InitArgs(xx,params);
02258 for (binz=0;binz<=nbinsz+1;binz++) {
02259 xx[2] = fZaxis.GetBinCenter(binz);
02260 for (biny=0;biny<=nbinsy+1;biny++) {
02261 xx[1] = fYaxis.GetBinCenter(biny);
02262 for (binx=0;binx<=nbinsx+1;binx++) {
02263 xx[0] = fXaxis.GetBinCenter(binx);
02264 if (!f1->IsInside(xx)) continue;
02265 TF1::RejectPoint(kFALSE);
02266 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
02267 Double_t error1 = GetBinError(bin);
02268 cu = c1*f1->EvalPar(xx);
02269 if (TF1::RejectedPoint()) continue;
02270 if (cu) w = GetBinContent(bin)/cu;
02271 else w = 0;
02272 SetBinContent(bin,w);
02273 if (fSumw2.fN) {
02274 if (cu != 0) fSumw2.fArray[bin] = error1*error1/(cu*cu);
02275 else fSumw2.fArray[bin] = 0;
02276 }
02277 }
02278 }
02279 }
02280 ResetStats();
02281 }
02282
02283
02284 void TH1::Divide(const TH1 *h1)
02285 {
02286
02287
02288
02289
02290
02291
02292
02293
02294
02295
02296
02297
02298
02299
02300
02301 if (!h1) {
02302 Error("Divide","Attempt to divide by a non-existing histogram");
02303 return;
02304 }
02305
02306 Int_t nbinsx = GetNbinsX();
02307 Int_t nbinsy = GetNbinsY();
02308 Int_t nbinsz = GetNbinsZ();
02309
02310
02311 try {
02312 CheckConsistency(this,h1);
02313 } catch(DifferentNumberOfBins&) {
02314 Error("Divide","Attempt to divide histograms with different number of bins");
02315 return;
02316 } catch(DifferentAxisLimits&) {
02317 Warning("Divide","Attempt to divide histograms with different axis limits");
02318 } catch(DifferentBinLimits&) {
02319 Warning("Divide","Attempt to divide histograms with different bin limits");
02320 }
02321
02322 if (fDimension < 2) nbinsy = -1;
02323 if (fDimension < 3) nbinsz = -1;
02324
02325
02326 if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
02327
02328
02329
02330
02331 ResetBit(kCanRebin);
02332
02333
02334 Int_t bin, binx, biny, binz;
02335 Double_t c0,c1,w;
02336 for (binz=0;binz<=nbinsz+1;binz++) {
02337 for (biny=0;biny<=nbinsy+1;biny++) {
02338 for (binx=0;binx<=nbinsx+1;binx++) {
02339 bin = GetBin(binx,biny,binz);
02340 c0 = GetBinContent(bin);
02341 c1 = h1->GetBinContent(bin);
02342 if (c1) w = c0/c1;
02343 else w = 0;
02344 SetBinContent(bin,w);
02345 if (fSumw2.fN) {
02346 Double_t e0 = GetBinError(bin);
02347 Double_t e1 = h1->GetBinError(bin);
02348 Double_t c12= c1*c1;
02349 if (!c1) { fSumw2.fArray[bin] = 0; continue;}
02350 fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0)/(c12*c12);
02351 }
02352 }
02353 }
02354 }
02355 ResetStats();
02356 }
02357
02358
02359
02360 void TH1::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
02361 {
02362
02363
02364
02365
02366
02367
02368
02369
02370
02371
02372
02373
02374
02375
02376
02377
02378
02379
02380
02381
02382
02383
02384 TString opt = option;
02385 opt.ToLower();
02386 Bool_t binomial = kFALSE;
02387 if (opt.Contains("b")) binomial = kTRUE;
02388 if (!h1 || !h2) {
02389 Error("Divide","Attempt to divide by a non-existing histogram");
02390 return;
02391 }
02392
02393 Int_t nbinsx = GetNbinsX();
02394 Int_t nbinsy = GetNbinsY();
02395 Int_t nbinsz = GetNbinsZ();
02396
02397 try {
02398 CheckConsistency(h1,h2);
02399 CheckConsistency(this,h1);
02400 } catch(DifferentNumberOfBins&) {
02401 Error("Divide","Attempt to divide histograms with different number of bins");
02402 return;
02403 } catch(DifferentAxisLimits&) {
02404 Warning("Divide","Attempt to divide histograms with different axis limits");
02405 } catch(DifferentBinLimits&) {
02406 Warning("Divide","Attempt to divide histograms with different bin limits");
02407 }
02408
02409 if (!c2) {
02410 Error("Divide","Coefficient of dividing histogram cannot be zero");
02411 return;
02412 }
02413
02414 if (fDimension < 2) nbinsy = -1;
02415 if (fDimension < 3) nbinsz = -1;
02416
02417
02418 if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
02419
02420 SetMinimum();
02421 SetMaximum();
02422
02423
02424
02425 ResetBit(kCanRebin);
02426
02427
02428 Int_t bin, binx, biny, binz;
02429 Double_t b1,b2,w,d1,d2;
02430 d1 = c1*c1;
02431 d2 = c2*c2;
02432 for (binz=0;binz<=nbinsz+1;binz++) {
02433 for (biny=0;biny<=nbinsy+1;biny++) {
02434 for (binx=0;binx<=nbinsx+1;binx++) {
02435 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
02436 b1 = h1->GetBinContent(bin);
02437 b2 = h2->GetBinContent(bin);
02438 if (b2) w = c1*b1/(c2*b2);
02439 else w = 0;
02440 SetBinContent(bin,w);
02441 if (fSumw2.fN) {
02442 Double_t e1 = h1->GetBinError(bin);
02443 Double_t e2 = h2->GetBinError(bin);
02444 Double_t b22= b2*b2*d2;
02445 if (!b2) { fSumw2.fArray[bin] = 0; continue;}
02446 if (binomial) {
02447 if (b1 != b2) {
02448
02449 w = b1/b2;
02450
02451
02452
02453 fSumw2.fArray[bin] = TMath::Abs( ( (1.-2.*w)*e1*e1 + w*w*e2*e2 )/(b2*b2) );
02454 } else {
02455
02456
02457 fSumw2.fArray[bin] = 0;
02458 }
02459 } else {
02460 fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22);
02461 }
02462 }
02463 }
02464 }
02465 }
02466 ResetStats();
02467 if (binomial)
02468
02469 SetEntries ( h2->GetEntries() );
02470 }
02471
02472
02473 void TH1::Draw(Option_t *option)
02474 {
02475
02476
02477
02478
02479
02480
02481
02482
02483
02484
02485
02486
02487
02488
02489
02490
02491
02492
02493
02494
02495
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508 TString opt = option;
02509 opt.ToLower();
02510 if (gPad) {
02511 if (!gPad->IsEditable()) gROOT->MakeDefCanvas();
02512 if (opt.Contains("same")) {
02513 if (gPad->GetX1() == 0 && gPad->GetX2() == 1 &&
02514 gPad->GetY1() == 0 && gPad->GetY2() == 1 &&
02515 gPad->GetListOfPrimitives()->GetSize()==0) opt.ReplaceAll("same","");
02516 } else {
02517
02518
02519 if (TestBit(kCanDelete)) gPad->GetListOfPrimitives()->Remove(this);
02520 gPad->Clear();
02521 }
02522 } else {
02523 if (opt.Contains("same")) opt.ReplaceAll("same","");
02524 }
02525 AppendPad(option);
02526 }
02527
02528
02529 TH1 *TH1::DrawCopy(Option_t *) const
02530 {
02531
02532
02533
02534
02535
02536
02537
02538
02539
02540
02541 AbstractMethod("DrawCopy");
02542 return 0;
02543 }
02544
02545
02546 TH1 *TH1::DrawNormalized(Option_t *option, Double_t norm) const
02547 {
02548
02549
02550
02551
02552
02553
02554
02555
02556
02557
02558
02559
02560
02561
02562
02563
02564
02565
02566 Double_t sum = GetSumOfWeights();
02567 if (sum == 0) {
02568 Error("DrawNormalized","Sum of weights is null. Cannot normalize histogram: %s",GetName());
02569 return 0;
02570 }
02571 Bool_t addStatus = TH1::AddDirectoryStatus();
02572 TH1::AddDirectory(kFALSE);
02573 TH1 *h = (TH1*)Clone();
02574 h->SetBit(kCanDelete);
02575 h->Scale(norm/sum);
02576 if (TMath::Abs(fMaximum+1111) > 1e-3) h->SetMaximum(fMaximum*norm/sum);
02577 if (TMath::Abs(fMinimum+1111) > 1e-3) h->SetMinimum(fMinimum*norm/sum);
02578 h->Draw(option);
02579 TH1::AddDirectory(addStatus);
02580 return h;
02581 }
02582
02583
02584 void TH1::DrawPanel()
02585 {
02586
02587
02588
02589
02590
02591 if (!fPainter) {Draw(); if (gPad) gPad->Update();}
02592 if (fPainter) fPainter->DrawPanel();
02593 }
02594
02595
02596 void TH1::Eval(TF1 *f1, Option_t *option)
02597 {
02598
02599
02600
02601
02602
02603
02604
02605
02606
02607
02608
02609
02610 Double_t x[3];
02611 Int_t range,stat,add,bin,binx,biny,binz,nbinsx, nbinsy, nbinsz;
02612 if (!f1) return;
02613 Double_t fu;
02614 Double_t e=0;
02615 TString opt = option;
02616 opt.ToLower();
02617 if (opt.Contains("a")) add = 1;
02618 else add = 0;
02619 if (opt.Contains("s")) stat = 1;
02620 else stat = 0;
02621 if (opt.Contains("r")) range = 1;
02622 else range = 0;
02623 nbinsx = fXaxis.GetNbins();
02624 nbinsy = fYaxis.GetNbins();
02625 nbinsz = fZaxis.GetNbins();
02626 if (!add) Reset();
02627
02628 for (binz=1;binz<=nbinsz;binz++) {
02629 x[2] = fZaxis.GetBinCenter(binz);
02630 for (biny=1;biny<=nbinsy;biny++) {
02631 x[1] = fYaxis.GetBinCenter(biny);
02632 for (binx=1;binx<=nbinsx;binx++) {
02633 bin = GetBin(binx,biny,binz);
02634 x[0] = fXaxis.GetBinCenter(binx);
02635 if (range && !f1->IsInside(x)) continue;
02636 fu = f1->Eval(x[0],x[1],x[2]);
02637 if (stat) fu = gRandom->PoissonD(fu);
02638 if (fSumw2.fN) e = fSumw2.fArray[bin];
02639 AddBinContent(bin,fu);
02640 if (fSumw2.fN) fSumw2.fArray[bin] = e+ TMath::Abs(fu);
02641 }
02642 }
02643 }
02644 }
02645
02646
02647 void TH1::ExecuteEvent(Int_t event, Int_t px, Int_t py)
02648 {
02649
02650
02651
02652
02653
02654
02655
02656
02657
02658 if (fPainter) fPainter->ExecuteEvent(event, px, py);
02659 }
02660
02661
02662 TH1* TH1::FFT(TH1* h_output, Option_t *option)
02663 {
02664
02665
02666
02667
02668
02669
02670
02671
02672
02673
02674
02675
02676
02677
02678
02679
02680
02681
02682
02683
02684
02685
02686
02687
02688
02689
02690
02691
02692
02693
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704
02705
02706 Int_t ndim[3];
02707 ndim[0] = this->GetNbinsX();
02708 ndim[1] = this->GetNbinsY();
02709 ndim[2] = this->GetNbinsZ();
02710
02711 TVirtualFFT *fft;
02712 TString opt = option;
02713 opt.ToUpper();
02714 if (!opt.Contains("2R")){
02715 if (!opt.Contains("2C") && !opt.Contains("2HC") && !opt.Contains("DHT")) {
02716
02717 opt.Append("R2C");
02718 }
02719 fft = TVirtualFFT::FFT(this->GetDimension(), ndim, opt.Data());
02720 }
02721 else {
02722
02723 Int_t ind = opt.Index("R2R", 3);
02724 Int_t *kind = new Int_t[2];
02725 char t;
02726 t = opt[ind+4];
02727 kind[0] = atoi(&t);
02728 if (h_output->GetDimension()>1) {
02729 t = opt[ind+5];
02730 kind[1] = atoi(&t);
02731 }
02732 fft = TVirtualFFT::SineCosine(this->GetDimension(), ndim, kind, option);
02733 delete [] kind;
02734 }
02735
02736 if (!fft) return 0;
02737 Int_t in=0;
02738 for (Int_t binx = 1; binx<=ndim[0]; binx++) {
02739 for (Int_t biny=1; biny<=ndim[1]; biny++) {
02740 for (Int_t binz=1; binz<=ndim[2]; binz++) {
02741 fft->SetPoint(in, this->GetBinContent(binx, biny, binz));
02742 in++;
02743 }
02744 }
02745 }
02746 fft->Transform();
02747 h_output = TransformHisto(fft, h_output, option);
02748 return h_output;
02749 }
02750
02751
02752 Int_t TH1::Fill(Double_t x)
02753 {
02754
02755
02756
02757
02758
02759
02760
02761
02762
02763
02764
02765
02766 if (fBuffer) return BufferFill(x,1);
02767
02768 Int_t bin;
02769 fEntries++;
02770 bin =fXaxis.FindBin(x);
02771 if (bin <0) return -1;
02772 AddBinContent(bin);
02773 if (fSumw2.fN) ++fSumw2.fArray[bin];
02774 if (bin == 0 || bin > fXaxis.GetNbins()) {
02775 if (!fgStatOverflows) return -1;
02776 }
02777 ++fTsumw;
02778 ++fTsumw2;
02779 fTsumwx += x;
02780 fTsumwx2 += x*x;
02781 return bin;
02782 }
02783
02784
02785 Int_t TH1::Fill(Double_t x, Double_t w)
02786 {
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797
02798
02799 if (fBuffer) return BufferFill(x,w);
02800
02801 Int_t bin;
02802 fEntries++;
02803 bin =fXaxis.FindBin(x);
02804 if (bin <0) return -1;
02805 AddBinContent(bin, w);
02806 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
02807 if (bin == 0 || bin > fXaxis.GetNbins()) {
02808 if (!fgStatOverflows) return -1;
02809 }
02810 Double_t z= (w > 0 ? w : -w);
02811 fTsumw += z;
02812 fTsumw2 += z*z;
02813 fTsumwx += z*x;
02814 fTsumwx2 += z*x*x;
02815 return bin;
02816 }
02817
02818
02819 Int_t TH1::Fill(const char *namex, Double_t w)
02820 {
02821
02822
02823
02824
02825
02826
02827
02828
02829
02830
02831 Int_t bin;
02832 fEntries++;
02833 bin =fXaxis.FindBin(namex);
02834 if (bin <0) return -1;
02835 AddBinContent(bin, w);
02836 if (fSumw2.fN) fSumw2.fArray[bin] += w*w;
02837 if (bin == 0 || bin > fXaxis.GetNbins()) return -1;
02838 Double_t x = fXaxis.GetBinCenter(bin);
02839 Double_t z= (w > 0 ? w : -w);
02840 fTsumw += z;
02841 fTsumw2 += z*z;
02842 fTsumwx += z*x;
02843 fTsumwx2 += z*x*x;
02844 return bin;
02845 }
02846
02847
02848 void TH1::FillN(Int_t ntimes, const Double_t *x, const Double_t *w, Int_t stride)
02849 {
02850
02851
02852
02853
02854
02855
02856
02857
02858
02859
02860
02861
02862
02863
02864
02865 Int_t bin,i;
02866
02867
02868
02869
02870
02871
02872
02873
02874
02875 fEntries += ntimes;
02876 Double_t ww = 1;
02877 Int_t nbins = fXaxis.GetNbins();
02878 ntimes *= stride;
02879 for (i=0;i<ntimes;i+=stride) {
02880 bin =fXaxis.FindBin(x[i]);
02881 if (bin <0) continue;
02882 if (w) ww = w[i];
02883 AddBinContent(bin, ww);
02884 if (fSumw2.fN) fSumw2.fArray[bin] += ww*ww;
02885 if (bin == 0 || bin > nbins) {
02886 if (!fgStatOverflows) continue;
02887 }
02888 Double_t z= (ww > 0 ? ww : -ww);
02889 fTsumw += z;
02890 fTsumw2 += z*z;
02891 fTsumwx += z*x[i];
02892 fTsumwx2 += z*x[i]*x[i];
02893 }
02894 }
02895
02896
02897 void TH1::FillRandom(const char *fname, Int_t ntimes)
02898 {
02899
02900
02901
02902
02903
02904
02905
02906
02907
02908
02909
02910
02911
02912
02913
02914
02915 Int_t bin, binx, ibin, loop;
02916 Double_t r1, x;
02917
02918 TF1 *f1 = (TF1*)gROOT->GetFunction(fname);
02919 if (!f1) { Error("FillRandom", "Unknown function: %s",fname); return; }
02920
02921
02922 Int_t first = fXaxis.GetFirst();
02923 Int_t last = fXaxis.GetLast();
02924 Int_t nbinsx = last-first+1;
02925
02926 Double_t *integral = new Double_t[nbinsx+1];
02927 integral[0] = 0;
02928 for (binx=1;binx<=nbinsx;binx++) {
02929 Double_t fint = f1->Integral(fXaxis.GetBinLowEdge(binx+first-1),fXaxis.GetBinUpEdge(binx+first-1));
02930 integral[binx] = integral[binx-1] + fint;
02931 }
02932
02933
02934 if (integral[nbinsx] == 0 ) {
02935 delete [] integral;
02936 Error("FillRandom", "Integral = zero"); return;
02937 }
02938 for (bin=1;bin<=nbinsx;bin++) integral[bin] /= integral[nbinsx];
02939
02940
02941 for (loop=0;loop<ntimes;loop++) {
02942 r1 = gRandom->Rndm(loop);
02943 ibin = TMath::BinarySearch(nbinsx,&integral[0],r1);
02944
02945
02946 x = fXaxis.GetBinLowEdge(ibin+first)
02947 +fXaxis.GetBinWidth(ibin+first)*(r1-integral[ibin])/(integral[ibin+1] - integral[ibin]);
02948 Fill(x, 1.);
02949 }
02950 delete [] integral;
02951 }
02952
02953
02954 void TH1::FillRandom(TH1 *h, Int_t ntimes)
02955 {
02956
02957
02958
02959
02960
02961
02962
02963
02964
02965
02966
02967
02968
02969
02970
02971
02972
02973
02974 if (!h) { Error("FillRandom", "Null histogram"); return; }
02975 if (fDimension != h->GetDimension()) {
02976 Error("FillRandom", "Histograms with different dimensions"); return;
02977 }
02978
02979
02980
02981 Int_t first = fXaxis.GetFirst();
02982 Int_t last = fXaxis.GetLast();
02983 Int_t nbins = last-first+1;
02984 if (CheckConsistency(this,h) && ntimes > 10*nbins) {
02985 Double_t sumw = h->Integral(first,last);
02986 if (sumw == 0) return;
02987 Double_t sumgen = 0;
02988 for (Int_t bin=first;bin<=last;bin++) {
02989 Double_t mean = h->GetBinContent(bin)*ntimes/sumw;
02990 Double_t cont = (Double_t)gRandom->Poisson(mean);
02991 sumgen += cont;
02992 AddBinContent(bin,cont);
02993 if (fSumw2.fN) fSumw2.fArray[bin] += cont;
02994 }
02995
02996
02997
02998
02999 Int_t i;
03000 if (sumgen < ntimes) {
03001
03002 for (i = Int_t(sumgen+0.5); i < ntimes; ++i)
03003 {
03004 Double_t x = h->GetRandom();
03005 Fill(x);
03006 }
03007 }
03008 else if (sumgen > ntimes) {
03009
03010 i = Int_t(sumgen+0.5);
03011 while( i > ntimes) {
03012 Double_t x = h->GetRandom();
03013 Int_t ibin = fXaxis.FindBin(x);
03014 Double_t y = GetBinContent(ibin);
03015
03016 if (y > 0) {
03017 SetBinContent(ibin, y-1.);
03018 i--;
03019 }
03020 }
03021 }
03022
03023 ResetStats();
03024 return;
03025 }
03026
03027 if (h->ComputeIntegral() ==0) return;
03028 Int_t loop;
03029 Double_t x;
03030 for (loop=0;loop<ntimes;loop++) {
03031 x = h->GetRandom();
03032 Fill(x);
03033 }
03034 }
03035
03036
03037
03038 Int_t TH1::FindBin(Double_t x, Double_t y, Double_t z)
03039 {
03040
03041
03042
03043
03044
03045
03046
03047
03048
03049
03050
03051 if (GetDimension() < 2) {
03052 return fXaxis.FindBin(x);
03053 }
03054 if (GetDimension() < 3) {
03055 Int_t nx = fXaxis.GetNbins()+2;
03056 Int_t binx = fXaxis.FindBin(x);
03057 Int_t biny = fYaxis.FindBin(y);
03058 return binx + nx*biny;
03059 }
03060 if (GetDimension() < 4) {
03061 Int_t nx = fXaxis.GetNbins()+2;
03062 Int_t ny = fYaxis.GetNbins()+2;
03063 Int_t binx = fXaxis.FindBin(x);
03064 Int_t biny = fYaxis.FindBin(y);
03065 Int_t binz = fZaxis.FindBin(z);
03066 return binx + nx*(biny +ny*binz);
03067 }
03068 return -1;
03069 }
03070
03071
03072 Int_t TH1::FindFixBin(Double_t x, Double_t y, Double_t z) const
03073 {
03074
03075
03076
03077
03078
03079
03080
03081
03082
03083
03084
03085 if (GetDimension() < 2) {
03086 return fXaxis.FindFixBin(x);
03087 }
03088 if (GetDimension() < 3) {
03089 Int_t nx = fXaxis.GetNbins()+2;
03090 Int_t binx = fXaxis.FindFixBin(x);
03091 Int_t biny = fYaxis.FindFixBin(y);
03092 return binx + nx*biny;
03093 }
03094 if (GetDimension() < 4) {
03095 Int_t nx = fXaxis.GetNbins()+2;
03096 Int_t ny = fYaxis.GetNbins()+2;
03097 Int_t binx = fXaxis.FindFixBin(x);
03098 Int_t biny = fYaxis.FindFixBin(y);
03099 Int_t binz = fZaxis.FindFixBin(z);
03100 return binx + nx*(biny +ny*binz);
03101 }
03102 return -1;
03103 }
03104
03105
03106 Int_t TH1::FindFirstBinAbove(Double_t threshold, Int_t axis) const
03107 {
03108
03109
03110
03111 if (axis != 1) {
03112 Warning("FindFirstBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
03113 axis = 1;
03114 }
03115 Int_t nbins = fXaxis.GetNbins();
03116 for (Int_t bin=1;bin<=nbins;bin++) {
03117 if (GetBinContent(bin) > threshold) return bin;
03118 }
03119 return -1;
03120 }
03121
03122
03123
03124 Int_t TH1::FindLastBinAbove(Double_t threshold, Int_t axis) const
03125 {
03126
03127
03128
03129 if (axis != 1) {
03130 Warning("FindLastBinAbove","Invalid axis number : %d, axis x assumed\n",axis);
03131 axis = 1;
03132 }
03133 Int_t nbins = fXaxis.GetNbins();
03134 for (Int_t bin=nbins;bin>=1;bin--) {
03135 if (GetBinContent(bin) > threshold) return bin;
03136 }
03137 return -1;
03138 }
03139
03140
03141 TObject *TH1::FindObject(const char *name) const
03142 {
03143
03144
03145 if (fFunctions) return fFunctions->FindObject(name);
03146 return 0;
03147 }
03148
03149
03150 TObject *TH1::FindObject(const TObject *obj) const
03151 {
03152
03153
03154 if (fFunctions) return fFunctions->FindObject(obj);
03155 return 0;
03156 }
03157
03158
03159 TFitResultPtr TH1::Fit(const char *fname ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
03160 {
03161
03162
03163
03164
03165
03166
03167
03168
03169
03170
03171
03172 char *linear;
03173 linear= (char*)strstr(fname, "++");
03174 TF1 *f1=0;
03175 TF2 *f2=0;
03176 TF3 *f3=0;
03177 Int_t ndim=GetDimension();
03178 if (linear){
03179 if (ndim<2){
03180 f1=new TF1(fname, fname, xxmin, xxmax);
03181 return Fit(f1,option,goption,xxmin,xxmax);
03182 }
03183 else if (ndim<3){
03184 f2=new TF2(fname, fname);
03185 return Fit(f2,option,goption,xxmin,xxmax);
03186 }
03187 else{
03188 f3=new TF3(fname, fname);
03189 return Fit(f3,option,goption,xxmin,xxmax);
03190 }
03191 }
03192
03193 else{
03194 f1 = (TF1*)gROOT->GetFunction(fname);
03195 if (!f1) { Printf("Unknown function: %s",fname); return -1; }
03196 return Fit(f1,option,goption,xxmin,xxmax);
03197 }
03198 }
03199
03200
03201 TFitResultPtr TH1::Fit(TF1 *f1 ,Option_t *option ,Option_t *goption, Double_t xxmin, Double_t xxmax)
03202 {
03203
03204
03205
03206
03207
03208
03209
03210
03211
03212
03213
03214
03215
03216
03217
03218
03219
03220
03221
03222
03223
03224
03225
03226
03227
03228
03229
03230
03231
03232
03233
03234
03235
03236
03237
03238
03239
03240
03241
03242
03243
03244
03245
03246
03247
03248
03249
03250
03251
03252
03253
03254
03255
03256
03257
03258
03259
03260
03261
03262
03263
03264
03265
03266
03267
03268
03269
03270
03271
03272
03273
03274
03275
03276
03277
03278
03279
03280
03281
03282
03283
03284
03285
03286
03287
03288
03289
03290
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303
03304
03305
03306
03307
03308
03309
03310
03311
03312
03313
03314
03315
03316
03317
03318
03319
03320
03321
03322
03323
03324
03325
03326
03327
03328
03329
03330
03331
03332
03333
03334
03335
03336
03337
03338
03339
03340
03341
03342
03343
03344
03345
03346
03347
03348
03349
03350
03351
03352
03353
03354
03355
03356
03357
03358
03359
03360
03361
03362
03363
03364
03365
03366
03367
03368
03369
03370
03371
03372
03373
03374
03375
03376
03377
03378
03379
03380
03381
03382
03383
03384
03385
03386
03387
03388
03389
03390
03391
03392
03393
03394
03395
03396
03397
03398
03399
03400
03401
03402
03403
03404
03405
03406
03407
03408
03409
03410
03411
03412
03413
03414
03415
03416
03417
03418
03419
03420
03421
03422
03423 Foption_t fitOption;
03424
03425 if (!FitOptionsMake(option,fitOption)) return 0;
03426
03427 ROOT::Fit::DataRange range(xxmin,xxmax);
03428 ROOT::Math::MinimizerOptions minOption;
03429
03430 return ROOT::Fit::FitObject(this, f1 , fitOption , minOption, goption, range);
03431 }
03432
03433
03434 void TH1::FitPanel()
03435 {
03436
03437
03438
03439
03440
03441 if (!gPad)
03442 gROOT->MakeDefCanvas();
03443
03444 if (!gPad) {
03445 Error("FitPanel", "Unable to create a default canvas");
03446 return;
03447 }
03448
03449
03450
03451 TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TFitEditor");
03452 if (handler && handler->LoadPlugin() != -1) {
03453 if (handler->ExecPlugin(2, gPad, this) == 0)
03454 Error("FitPanel", "Unable to create the FitPanel");
03455 }
03456 else
03457 Error("FitPanel", "Unable to find the FitPanel plug-in");
03458 }
03459
03460
03461 TH1 *TH1::GetAsymmetry(TH1* h2, Double_t c2, Double_t dc2)
03462 {
03463
03464
03465
03466
03467
03468
03469
03470
03471
03472
03473
03474
03475
03476
03477
03478
03479
03480
03481
03482
03483
03484
03485
03486
03487
03488
03489
03490 Bool_t addStatus = TH1::AddDirectoryStatus();
03491 TH1::AddDirectory(kFALSE);
03492 TH1 *asym = (TH1*)Clone();
03493 asym->Sumw2();
03494 TH1 *top = (TH1*)asym->Clone();
03495 TH1 *bottom = (TH1*)asym->Clone();
03496 TH1::AddDirectory(addStatus);
03497
03498
03499 TH1 *h1 = this;
03500 top->Add(h1,h2,1,-c2);
03501 bottom->Add(h1,h2,1,c2);
03502 asym->Divide(top,bottom);
03503
03504 Int_t xmax = asym->GetNbinsX();
03505 Int_t ymax = asym->GetNbinsY();
03506 Int_t zmax = asym->GetNbinsZ();
03507 Double_t bot, error, a, b, da, db;
03508
03509
03510
03511 for(Int_t i=1; i<= xmax; i++){
03512 for(Int_t j=1; j<= ymax; j++){
03513 for(Int_t k=1; k<= zmax; k++){
03514
03515
03516
03517 a = h1->GetBinContent(i,j,k);
03518 b = h2->GetBinContent(i,j,k);
03519 bot = bottom->GetBinContent(i,j,k);
03520
03521
03522
03523
03524 if(bot < 1e-6){}
03525 else{
03526
03527 da = h1->GetBinError(i,j,k);
03528 db = h2->GetBinError(i,j,k);
03529 error = 2*TMath::Sqrt(a*a*c2*c2*db*db + c2*c2*b*b*da*da+a*a*b*b*dc2*dc2)/(bot*bot);
03530 asym->SetBinError(i,j,k,error);
03531 }
03532 }
03533 }
03534 }
03535 delete top;
03536 delete bottom;
03537 return asym;
03538 }
03539
03540
03541 Int_t TH1::GetDefaultBufferSize()
03542 {
03543
03544
03545
03546
03547 return fgBufferSize;
03548 }
03549
03550
03551
03552 Bool_t TH1::GetDefaultSumw2()
03553 {
03554
03555
03556
03557
03558 return fgDefaultSumw2;
03559 }
03560
03561
03562
03563 Double_t TH1::GetEntries() const
03564 {
03565
03566
03567 if (fBuffer) ((TH1*)this)->BufferEmpty();
03568
03569 return fEntries;
03570 }
03571
03572
03573 Double_t TH1::GetEffectiveEntries() const
03574 {
03575
03576
03577
03578
03579
03580 Stat_t s[kNstat];
03581 this->GetStats(s);
03582 return (s[1] ? s[0]*s[0]/s[1] : TMath::Abs(s[0]) );
03583 }
03584
03585
03586 char *TH1::GetObjectInfo(Int_t px, Int_t py) const
03587 {
03588
03589
03590
03591
03592 return ((TH1*)this)->GetPainter()->GetObjectInfo(px,py);
03593 }
03594
03595
03596 TVirtualHistPainter *TH1::GetPainter(Option_t *option)
03597 {
03598
03599
03600 if (!fPainter) {
03601 TString opt = option;
03602 opt.ToLower();
03603 if (opt.Contains("gl") || gStyle->GetCanvasPreferGL()) {
03604
03605 TPluginHandler *handler = gROOT->GetPluginManager()->FindHandler("TGLHistPainter");
03606
03607 if (handler && handler->LoadPlugin() != -1)
03608 fPainter = reinterpret_cast<TVirtualHistPainter *>(handler->ExecPlugin(1, this));
03609 }
03610 }
03611
03612 if (!fPainter) fPainter = TVirtualHistPainter::HistPainter(this);
03613
03614 return fPainter;
03615 }
03616
03617
03618 Int_t TH1::GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum)
03619 {
03620
03621
03622
03623
03624
03625
03626
03627
03628
03629
03630
03631
03632
03633
03634
03635
03636
03637
03638
03639
03640
03641
03642
03643
03644
03645
03646
03647
03648
03649
03650
03651
03652
03653
03654
03655
03656
03657
03658
03659
03660
03661
03662
03663
03664
03665
03666
03667
03668
03669
03670
03671
03672
03673
03674
03675
03676
03677
03678
03679
03680
03681
03682
03683
03684
03685 if (GetDimension() > 1) {
03686 Error("GetQuantiles","Only available for 1-d histograms");
03687 return 0;
03688 }
03689
03690 const Int_t nbins = GetXaxis()->GetNbins();
03691 if (!fIntegral) ComputeIntegral();
03692 if (fIntegral[nbins+1] != fEntries) ComputeIntegral();
03693
03694 Int_t i, ibin;
03695 Double_t *prob = (Double_t*)probSum;
03696 Int_t nq = nprobSum;
03697 if (probSum == 0) {
03698 nq = nbins+1;
03699 prob = new Double_t[nq];
03700 prob[0] = 0;
03701 for (i=1;i<nq;i++) {
03702 prob[i] = fIntegral[i]/fIntegral[nbins];
03703 }
03704 }
03705
03706 for (i = 0; i < nq; i++) {
03707 ibin = TMath::BinarySearch(nbins,fIntegral,prob[i]);
03708 while (ibin < nbins-1 && fIntegral[ibin+1] == prob[i]) {
03709 if (fIntegral[ibin+2] == prob[i]) ibin++;
03710 else break;
03711 }
03712 q[i] = GetBinLowEdge(ibin+1);
03713 const Double_t dint = fIntegral[ibin+1]-fIntegral[ibin];
03714 if (dint > 0) q[i] += GetBinWidth(ibin+1)*(prob[i]-fIntegral[ibin])/dint;
03715 }
03716
03717 if (!probSum) delete [] prob;
03718 return nq;
03719 }
03720
03721
03722 Int_t TH1::FitOptionsMake(Option_t *choptin, Foption_t &fitOption)
03723 {
03724
03725
03726
03727 if (choptin == 0) return 1;
03728 if (strlen(choptin) == 0) return 1;
03729 TString opt = choptin;
03730 opt.ToUpper();
03731 if (opt.Contains("Q")) fitOption.Quiet = 1;
03732 if (opt.Contains("V")) {fitOption.Verbose = 1; fitOption.Quiet = 0;}
03733 if (opt.Contains("X")) fitOption.Chi2 = 1;
03734 if (opt.Contains("L")) fitOption.Like = 1;
03735 if (opt.Contains("LL")) fitOption.Like = 2;
03736 if (opt.Contains("W")) fitOption.W1 = 1;
03737 if (opt.Contains("WW")) fitOption.W1 = 2;
03738 if (opt.Contains("E")) fitOption.Errors = 1;
03739 if (opt.Contains("M")) fitOption.More = 1;
03740 if (opt.Contains("R")) fitOption.Range = 1;
03741 if (opt.Contains("G")) fitOption.Gradient= 1;
03742 if (opt.Contains("N")) fitOption.Nostore = 1;
03743 if (opt.Contains("0")) fitOption.Nograph = 1;
03744 if (opt.Contains("+")) fitOption.Plus = 1;
03745 if (opt.Contains("I")) fitOption.Integral= 1;
03746 if (opt.Contains("B")) fitOption.Bound = 1;
03747 if (opt.Contains("U")) {fitOption.User = 1; fitOption.Like = 0;}
03748 if (opt.Contains("F")) fitOption.Minuit = 1;
03749 if (opt.Contains("C")) fitOption.Nochisq = 1;
03750 if (opt.Contains("S")) fitOption.StoreResult = 1;
03751 return 1;
03752 }
03753
03754
03755 void H1InitGaus()
03756 {
03757
03758
03759
03760 Double_t allcha, sumx, sumx2, x, val, rms, mean;
03761 Int_t bin;
03762 const Double_t sqrtpi = 2.506628;
03763
03764
03765 TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
03766 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
03767 Int_t hxfirst = hFitter->GetXfirst();
03768 Int_t hxlast = hFitter->GetXlast();
03769 Double_t valmax = curHist->GetBinContent(hxfirst);
03770 Double_t binwidx = curHist->GetBinWidth(hxfirst);
03771 allcha = sumx = sumx2 = 0;
03772 for (bin=hxfirst;bin<=hxlast;bin++) {
03773 x = curHist->GetBinCenter(bin);
03774 val = TMath::Abs(curHist->GetBinContent(bin));
03775 if (val > valmax) valmax = val;
03776 sumx += val*x;
03777 sumx2 += val*x*x;
03778 allcha += val;
03779 }
03780 if (allcha == 0) return;
03781 mean = sumx/allcha;
03782 rms = sumx2/allcha - mean*mean;
03783 if (rms > 0) rms = TMath::Sqrt(rms);
03784 else rms = 0;
03785 if (rms == 0) rms = binwidx*(hxlast-hxfirst+1)/4;
03786
03787
03788
03789
03790
03791
03792 Double_t constant = 0.5*(valmax+binwidx*allcha/(sqrtpi*rms));
03793
03794
03795
03796
03797
03798 Double_t xmin = curHist->GetXaxis()->GetXmin();
03799 Double_t xmax = curHist->GetXaxis()->GetXmax();
03800 if ((mean < xmin || mean > xmax) && rms > (xmax-xmin)) {
03801 mean = 0.5*(xmax+xmin);
03802 rms = 0.5*(xmax-xmin);
03803 }
03804 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
03805 f1->SetParameter(0,constant);
03806 f1->SetParameter(1,mean);
03807 f1->SetParameter(2,rms);
03808 f1->SetParLimits(2,0,10*rms);
03809 }
03810
03811
03812 void H1InitExpo()
03813 {
03814
03815
03816
03817 Double_t constant, slope;
03818 Int_t ifail;
03819 TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
03820 Int_t hxfirst = hFitter->GetXfirst();
03821 Int_t hxlast = hFitter->GetXlast();
03822 Int_t nchanx = hxlast - hxfirst + 1;
03823
03824 H1LeastSquareLinearFit(-nchanx, constant, slope, ifail);
03825
03826 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
03827 f1->SetParameter(0,constant);
03828 f1->SetParameter(1,slope);
03829
03830 }
03831
03832
03833 void H1InitPolynom()
03834 {
03835
03836
03837
03838 Double_t fitpar[25];
03839
03840 TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
03841 TF1 *f1 = (TF1*)hFitter->GetUserFunc();
03842 Int_t hxfirst = hFitter->GetXfirst();
03843 Int_t hxlast = hFitter->GetXlast();
03844 Int_t nchanx = hxlast - hxfirst + 1;
03845 Int_t npar = f1->GetNpar();
03846
03847 if (nchanx <=1 || npar == 1) {
03848 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
03849 fitpar[0] = curHist->GetSumOfWeights()/Double_t(nchanx);
03850 } else {
03851 H1LeastSquareFit( nchanx, npar, fitpar);
03852 }
03853 for (Int_t i=0;i<npar;i++) f1->SetParameter(i, fitpar[i]);
03854 }
03855
03856
03857 void H1LeastSquareFit(Int_t n, Int_t m, Double_t *a)
03858 {
03859
03860
03861
03862
03863
03864
03865
03866
03867
03868
03869
03870 const Double_t zero = 0.;
03871 const Double_t one = 1.;
03872 const Int_t idim = 20;
03873
03874 Double_t b[400] ;
03875 Int_t i, k, l, ifail;
03876 Double_t power;
03877 Double_t da[20], xk, yk;
03878
03879 if (m <= 2) {
03880 H1LeastSquareLinearFit(n, a[0], a[1], ifail);
03881 return;
03882 }
03883 if (m > idim || m > n) return;
03884 b[0] = Double_t(n);
03885 da[0] = zero;
03886 for (l = 2; l <= m; ++l) {
03887 b[l-1] = zero;
03888 b[m + l*20 - 21] = zero;
03889 da[l-1] = zero;
03890 }
03891 TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
03892 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
03893 Int_t hxfirst = hFitter->GetXfirst();
03894 Int_t hxlast = hFitter->GetXlast();
03895 for (k = hxfirst; k <= hxlast; ++k) {
03896 xk = curHist->GetBinCenter(k);
03897 yk = curHist->GetBinContent(k);
03898 power = one;
03899 da[0] += yk;
03900 for (l = 2; l <= m; ++l) {
03901 power *= xk;
03902 b[l-1] += power;
03903 da[l-1] += power*yk;
03904 }
03905 for (l = 2; l <= m; ++l) {
03906 power *= xk;
03907 b[m + l*20 - 21] += power;
03908 }
03909 }
03910 for (i = 3; i <= m; ++i) {
03911 for (k = i; k <= m; ++k) {
03912 b[k - 1 + (i-1)*20 - 21] = b[k + (i-2)*20 - 21];
03913 }
03914 }
03915 H1LeastSquareSeqnd(m, b, idim, ifail, 1, da);
03916
03917 for (i=0; i<m; ++i) a[i] = da[i];
03918
03919 }
03920
03921
03922 void H1LeastSquareLinearFit(Int_t ndata, Double_t &a0, Double_t &a1, Int_t &ifail)
03923 {
03924
03925
03926
03927
03928
03929
03930
03931
03932 Double_t xbar, ybar, x2bar;
03933 Int_t i, n;
03934 Double_t xybar;
03935 Double_t fn, xk, yk;
03936 Double_t det;
03937
03938 n = TMath::Abs(ndata);
03939 ifail = -2;
03940 xbar = ybar = x2bar = xybar = 0;
03941 TVirtualFitter *hFitter = TVirtualFitter::GetFitter();
03942 TH1 *curHist = (TH1*)hFitter->GetObjectFit();
03943 Int_t hxfirst = hFitter->GetXfirst();
03944 Int_t hxlast = hFitter->GetXlast();
03945 for (i = hxfirst; i <= hxlast; ++i) {
03946 xk = curHist->GetBinCenter(i);
03947 yk = curHist->GetBinContent(i);
03948 if (ndata < 0) {
03949 if (yk <= 0) yk = 1e-9;
03950 yk = TMath::Log(yk);
03951 }
03952 xbar += xk;
03953 ybar += yk;
03954 x2bar += xk*xk;
03955 xybar += xk*yk;
03956 }
03957 fn = Double_t(n);
03958 det = fn*x2bar - xbar*xbar;
03959 ifail = -1;
03960 if (det <= 0) {
03961 a0 = ybar/fn;
03962 a1 = 0;
03963 return;
03964 }
03965 ifail = 0;
03966 a0 = (x2bar*ybar - xbar*xybar) / det;
03967 a1 = (fn*xybar - xbar*ybar) / det;
03968
03969 }
03970
03971
03972 void H1LeastSquareSeqnd(Int_t n, Double_t *a, Int_t idim, Int_t &ifail, Int_t k, Double_t *b)
03973 {
03974
03975
03976
03977
03978
03979
03980 Int_t a_dim1, a_offset, b_dim1, b_offset;
03981 Int_t nmjp1, i, j, l;
03982 Int_t im1, jp1, nm1, nmi;
03983 Double_t s1, s21, s22;
03984 const Double_t one = 1.;
03985
03986
03987 b_dim1 = idim;
03988 b_offset = b_dim1 + 1;
03989 b -= b_offset;
03990 a_dim1 = idim;
03991 a_offset = a_dim1 + 1;
03992 a -= a_offset;
03993
03994 if (idim < n) return;
03995
03996 ifail = 0;
03997 for (j = 1; j <= n; ++j) {
03998 if (a[j + j*a_dim1] <= 0) { ifail = -1; return; }
03999 a[j + j*a_dim1] = one / a[j + j*a_dim1];
04000 if (j == n) continue;
04001 jp1 = j + 1;
04002 for (l = jp1; l <= n; ++l) {
04003 a[j + l*a_dim1] = a[j + j*a_dim1] * a[l + j*a_dim1];
04004 s1 = -a[l + (j+1)*a_dim1];
04005 for (i = 1; i <= j; ++i) { s1 = a[l + i*a_dim1] * a[i + (j+1)*a_dim1] + s1; }
04006 a[l + (j+1)*a_dim1] = -s1;
04007 }
04008 }
04009 if (k <= 0) return;
04010
04011 for (l = 1; l <= k; ++l) {
04012 b[l*b_dim1 + 1] = a[a_dim1 + 1]*b[l*b_dim1 + 1];
04013 }
04014 if (n == 1) return;
04015 for (l = 1; l <= k; ++l) {
04016 for (i = 2; i <= n; ++i) {
04017 im1 = i - 1;
04018 s21 = -b[i + l*b_dim1];
04019 for (j = 1; j <= im1; ++j) {
04020 s21 = a[i + j*a_dim1]*b[j + l*b_dim1] + s21;
04021 }
04022 b[i + l*b_dim1] = -a[i + i*a_dim1]*s21;
04023 }
04024 nm1 = n - 1;
04025 for (i = 1; i <= nm1; ++i) {
04026 nmi = n - i;
04027 s22 = -b[nmi + l*b_dim1];
04028 for (j = 1; j <= i; ++j) {
04029 nmjp1 = n - j + 1;
04030 s22 = a[nmi + nmjp1*a_dim1]*b[nmjp1 + l*b_dim1] + s22;
04031 }
04032 b[nmi + l*b_dim1] = -s22;
04033 }
04034 }
04035 }
04036
04037
04038
04039 Int_t TH1::GetBin(Int_t binx, Int_t biny, Int_t binz) const
04040 {
04041
04042
04043
04044
04045
04046
04047
04048
04049
04050
04051
04052
04053
04054
04055
04056
04057
04058
04059
04060
04061
04062
04063
04064
04065 Int_t nx, ny, nz;
04066 if (GetDimension() < 2) {
04067 nx = fXaxis.GetNbins()+2;
04068 if (binx < 0) binx = 0;
04069 if (binx >= nx) binx = nx-1;
04070 return binx;
04071 }
04072 if (GetDimension() < 3) {
04073 nx = fXaxis.GetNbins()+2;
04074 if (binx < 0) binx = 0;
04075 if (binx >= nx) binx = nx-1;
04076 ny = fYaxis.GetNbins()+2;
04077 if (biny < 0) biny = 0;
04078 if (biny >= ny) biny = ny-1;
04079 return binx + nx*biny;
04080 }
04081 if (GetDimension() < 4) {
04082 nx = fXaxis.GetNbins()+2;
04083 if (binx < 0) binx = 0;
04084 if (binx >= nx) binx = nx-1;
04085 ny = fYaxis.GetNbins()+2;
04086 if (biny < 0) biny = 0;
04087 if (biny >= ny) biny = ny-1;
04088 nz = fZaxis.GetNbins()+2;
04089 if (binz < 0) binz = 0;
04090 if (binz >= nz) binz = nz-1;
04091 return binx + nx*(biny +ny*binz);
04092 }
04093 return -1;
04094 }
04095
04096 void TH1::GetBinXYZ(Int_t binglobal, Int_t &binx, Int_t &biny, Int_t &binz) const
04097 {
04098
04099
04100
04101 Int_t nx = fXaxis.GetNbins()+2;
04102 Int_t ny = fYaxis.GetNbins()+2;
04103
04104 if (GetDimension() < 2) {
04105 binx = binglobal%nx;
04106 }
04107 if (GetDimension() < 3) {
04108 binx = binglobal%nx;
04109 biny = ((binglobal-binx)/nx)%ny;
04110 }
04111 if (GetDimension() < 4) {
04112 binx = binglobal%nx;
04113 biny = ((binglobal-binx)/nx)%ny;
04114 binz = ((binglobal-binx)/nx -biny)/ny;
04115 }
04116 }
04117
04118
04119
04120 Double_t TH1::GetRandom() const
04121 {
04122
04123
04124
04125
04126
04127
04128
04129 if (fDimension > 1) {
04130 Error("GetRandom","Function only valid for 1-d histograms");
04131 return 0;
04132 }
04133 Int_t nbinsx = GetNbinsX();
04134 Double_t integral;
04135 if (fIntegral) {
04136 if (fIntegral[nbinsx+1] != fEntries) integral = ((TH1*)this)->ComputeIntegral();
04137 } else {
04138 integral = ((TH1*)this)->ComputeIntegral();
04139 if (integral == 0 || fIntegral == 0) return 0;
04140 }
04141 Double_t r1 = gRandom->Rndm();
04142 Int_t ibin = TMath::BinarySearch(nbinsx,fIntegral,r1);
04143 Double_t x = GetBinLowEdge(ibin+1);
04144 if (r1 > fIntegral[ibin]) x +=
04145 GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1] - fIntegral[ibin]);
04146 return x;
04147 }
04148
04149
04150 Double_t TH1::GetBinContent(Int_t) const
04151 {
04152
04153
04154
04155
04156
04157
04158
04159
04160
04161
04162
04163
04164
04165
04166
04167
04168
04169 AbstractMethod("GetBinContent");
04170 return 0;
04171 }
04172
04173
04174 Double_t TH1::GetBinContent(Int_t binx, Int_t biny) const
04175 {
04176
04177
04178
04179
04180
04181 Int_t bin = GetBin(binx,biny);
04182 return GetBinContent(bin);
04183 }
04184
04185
04186 Double_t TH1::GetBinContent(Int_t binx, Int_t biny, Int_t binz) const
04187 {
04188
04189
04190
04191
04192
04193 Int_t bin = GetBin(binx,biny,binz);
04194 return GetBinContent(bin);
04195 }
04196
04197
04198 Double_t TH1::GetBinWithContent(Double_t c, Int_t &binx, Int_t firstx, Int_t lastx,Double_t maxdiff) const
04199 {
04200
04201
04202
04203
04204
04205
04206
04207
04208
04209
04210
04211
04212
04213 if (fDimension > 1) {
04214 binx = 0;
04215 Error("GetBinWithContent","function is only valid for 1-D histograms");
04216 return 0;
04217 }
04218 if (firstx <= 0) firstx = 1;
04219 if (lastx < firstx) lastx = fXaxis.GetNbins();
04220 Int_t binminx = 0;
04221 Double_t diff, curmax = 1.e240;
04222 for (Int_t i=firstx;i<=lastx;i++) {
04223 diff = TMath::Abs(GetBinContent(i)-c);
04224 if (diff <= 0) {binx = i; return diff;}
04225 if (diff < curmax && diff <= maxdiff) {curmax = diff, binminx=i;}
04226 }
04227 binx = binminx;
04228 return curmax;
04229 }
04230
04231
04232 TAxis *TH1::GetXaxis() const
04233 {
04234
04235
04236 return &((TH1*)this)->fXaxis;
04237 }
04238
04239
04240
04241 TAxis *TH1::GetYaxis() const
04242 {
04243
04244
04245 return &((TH1*)this)->fYaxis;
04246 }
04247
04248 TAxis *TH1::GetZaxis() const
04249 {
04250
04251
04252 return &((TH1*)this)->fZaxis;
04253 }
04254
04255
04256 Double_t TH1::Interpolate(Double_t x)
04257 {
04258
04259
04260
04261
04262 Int_t xbin = FindBin(x);
04263 Double_t x0,x1,y0,y1;
04264
04265 if(x<=GetBinCenter(1)) {
04266 return GetBinContent(1);
04267 } else if(x>=GetBinCenter(GetNbinsX())) {
04268 return GetBinContent(GetNbinsX());
04269 } else {
04270 if(x<=GetBinCenter(xbin)) {
04271 y0 = GetBinContent(xbin-1);
04272 x0 = GetBinCenter(xbin-1);
04273 y1 = GetBinContent(xbin);
04274 x1 = GetBinCenter(xbin);
04275 } else {
04276 y0 = GetBinContent(xbin);
04277 x0 = GetBinCenter(xbin);
04278 y1 = GetBinContent(xbin+1);
04279 x1 = GetBinCenter(xbin+1);
04280 }
04281 return y0 + (x-x0)*((y1-y0)/(x1-x0));
04282 }
04283 }
04284
04285
04286 Double_t TH1::Interpolate(Double_t, Double_t)
04287 {
04288
04289
04290 Error("Interpolate","This function must be called with 1 argument for a TH1");
04291 return 0;
04292 }
04293
04294
04295 Double_t TH1::Interpolate(Double_t, Double_t, Double_t)
04296 {
04297
04298
04299 Error("Interpolate","This function must be called with 1 argument for a TH1");
04300 return 0;
04301 }
04302
04303
04304 Bool_t TH1::IsBinOverflow(Int_t bin) const
04305 {
04306
04307
04308 Int_t binx, biny, binz;
04309 GetBinXYZ(bin, binx, biny, binz);
04310
04311 if ( fDimension == 1 )
04312 return binx >= GetNbinsX() + 1;
04313 else if ( fDimension == 2 )
04314 return (binx >= GetNbinsX() + 1) ||
04315 (biny >= GetNbinsY() + 1);
04316 else if ( fDimension == 3 )
04317 return (binx >= GetNbinsX() + 1) ||
04318 (biny >= GetNbinsY() + 1) ||
04319 (binz >= GetNbinsZ() + 1);
04320 else
04321 return 0;
04322 }
04323
04324
04325
04326 Bool_t TH1::IsBinUnderflow(Int_t bin) const
04327 {
04328
04329
04330 Int_t binx, biny, binz;
04331 GetBinXYZ(bin, binx, biny, binz);
04332
04333 if ( fDimension == 1 )
04334 return (binx <= 0);
04335 else if ( fDimension == 2 )
04336 return (binx <= 0 || biny <= 0);
04337 else if ( fDimension == 3 )
04338 return (binx <= 0 || biny <= 0 || binz <= 0);
04339 else
04340 return 0;
04341 }
04342
04343
04344 void TH1::LabelsDeflate(Option_t *ax)
04345 {
04346
04347
04348 Int_t iaxis = AxisChoice(ax);
04349 TAxis *axis = 0;
04350 if (iaxis == 1) axis = GetXaxis();
04351 if (iaxis == 2) axis = GetYaxis();
04352 if (iaxis == 3) axis = GetZaxis();
04353 if (!axis) return;
04354 if (!axis->GetLabels()) return;
04355 TIter next(axis->GetLabels());
04356 TObject *obj;
04357 Int_t nbins = 0;
04358 while ((obj = next())) {
04359 if (obj->GetUniqueID()) nbins++;
04360 }
04361 if (nbins < 1) nbins = 1;
04362 TH1 *hold = (TH1*)Clone();
04363 hold->SetDirectory(0);
04364
04365 Bool_t timedisp = axis->GetTimeDisplay();
04366 Double_t xmin = axis->GetXmin();
04367 Double_t xmax = axis->GetBinUpEdge(nbins);
04368 if (xmax <= xmin) xmax = xmin +nbins;
04369 axis->SetRange(0,0);
04370 axis->Set(nbins,xmin,xmax);
04371
04372
04373
04374 Int_t nbinsx = hold->GetXaxis()->GetNbins();
04375 Int_t nbinsy = hold->GetYaxis()->GetNbins();
04376 Int_t nbinsz = hold->GetZaxis()->GetNbins();
04377 Int_t ncells = nbinsx+2;
04378 if (GetDimension() > 1) ncells *= nbinsy+2;
04379 if (GetDimension() > 2) ncells *= nbinsz+2;
04380 SetBinsLength(ncells);
04381 Int_t errors = fSumw2.fN;
04382 if (errors) fSumw2.Set(ncells);
04383 axis->SetTimeDisplay(timedisp);
04384
04385
04386 Double_t err,cu;
04387 Int_t bin,ibin,binx,biny,binz;
04388 Double_t oldEntries = fEntries;
04389 for (binz=1;binz<=nbinsz;binz++) {
04390 for (biny=1;biny<=nbinsy;biny++) {
04391 for (binx=1;binx<=nbinsx;binx++) {
04392 bin = hold->GetBin(binx,biny,binz);
04393 ibin= GetBin(binx,biny,binz);
04394 cu = hold->GetBinContent(bin);
04395 SetBinContent(ibin,cu);
04396 if (errors) {
04397 err = hold->GetBinError(bin);
04398 SetBinError(ibin,err);
04399 }
04400 }
04401 }
04402 }
04403 fEntries = oldEntries;
04404 delete hold;
04405 }
04406
04407
04408 void TH1::LabelsInflate(Option_t *ax)
04409 {
04410
04411
04412
04413
04414 Int_t iaxis = AxisChoice(ax);
04415 TAxis *axis = 0;
04416 if (iaxis == 1) axis = GetXaxis();
04417 if (iaxis == 2) axis = GetYaxis();
04418 if (iaxis == 3) axis = GetZaxis();
04419 if (!axis) return;
04420
04421 TH1 *hold = (TH1*)Clone();
04422 hold->SetDirectory(0);
04423
04424 Bool_t timedisp = axis->GetTimeDisplay();
04425 Int_t nbxold = fXaxis.GetNbins();
04426 Int_t nbyold = fYaxis.GetNbins();
04427 Int_t nbzold = fZaxis.GetNbins();
04428 Int_t nbins = axis->GetNbins();
04429 Double_t xmin = axis->GetXmin();
04430 Double_t xmax = axis->GetXmax();
04431 xmax = xmin + 2*(xmax-xmin);
04432 axis->SetRange(0,0);
04433 axis->Set(2*nbins,xmin,xmax);
04434 Int_t nbinsx = fXaxis.GetNbins();
04435 Int_t nbinsy = fYaxis.GetNbins();
04436 Int_t nbinsz = fZaxis.GetNbins();
04437 Int_t ncells = nbinsx+2;
04438 if (GetDimension() > 1) ncells *= nbinsy+2;
04439 if (GetDimension() > 2) ncells *= nbinsz+2;
04440 SetBinsLength(ncells);
04441 Int_t errors = fSumw2.fN;
04442 if (errors) fSumw2.Set(ncells);
04443 axis->SetTimeDisplay(timedisp);
04444
04445
04446 Double_t err,cu;
04447 Double_t oldEntries = fEntries;
04448 Int_t bin,ibin,binx,biny,binz;
04449 Int_t binxmin = 1;
04450 Int_t binxmax = nbinsx;
04451 if (fDimension > 1) {binxmin--; binxmax++;}
04452 for (binz=1;binz<=nbinsz;binz++) {
04453 for (biny=1;biny<=nbinsy;biny++) {
04454 for (binx=binxmin;binx<=binxmax;binx++) {
04455 bin = hold->GetBin(binx,biny,binz);
04456 ibin= GetBin(binx,biny,binz);
04457 if (binx > nbxold+1 || biny > nbyold || binz > nbzold) bin = -1;
04458 if (bin > 0) cu = hold->GetBinContent(bin);
04459 else cu = 0;
04460 SetBinContent(ibin,cu);
04461 if (errors) {
04462 if (bin > 0) err = hold->GetBinError(bin);
04463 else err = 0;
04464 SetBinError(ibin,err);
04465 }
04466 }
04467 }
04468 }
04469 fEntries = oldEntries;
04470 delete hold;
04471 }
04472
04473
04474 void TH1::LabelsOption(Option_t *option, Option_t *ax)
04475 {
04476
04477
04478
04479
04480
04481
04482
04483
04484
04485 Int_t iaxis = AxisChoice(ax);
04486 TAxis *axis = 0;
04487 if (iaxis == 1) axis = GetXaxis();
04488 if (iaxis == 2) axis = GetYaxis();
04489 if (iaxis == 3) axis = GetZaxis();
04490 if (!axis) return;
04491 THashList *labels = axis->GetLabels();
04492 if (!labels) {
04493 Warning("LabelsOption","Cannot sort. No labels");
04494 return;
04495 }
04496 TString opt = option;
04497 opt.ToLower();
04498 if (opt.Contains("h")) {
04499 axis->SetBit(TAxis::kLabelsHori);
04500 axis->ResetBit(TAxis::kLabelsVert);
04501 axis->ResetBit(TAxis::kLabelsDown);
04502 axis->ResetBit(TAxis::kLabelsUp);
04503 }
04504 if (opt.Contains("v")) {
04505 axis->SetBit(TAxis::kLabelsVert);
04506 axis->ResetBit(TAxis::kLabelsHori);
04507 axis->ResetBit(TAxis::kLabelsDown);
04508 axis->ResetBit(TAxis::kLabelsUp);
04509 }
04510 if (opt.Contains("u")) {
04511 axis->SetBit(TAxis::kLabelsUp);
04512 axis->ResetBit(TAxis::kLabelsVert);
04513 axis->ResetBit(TAxis::kLabelsDown);
04514 axis->ResetBit(TAxis::kLabelsHori);
04515 }
04516 if (opt.Contains("d")) {
04517 axis->SetBit(TAxis::kLabelsDown);
04518 axis->ResetBit(TAxis::kLabelsVert);
04519 axis->ResetBit(TAxis::kLabelsHori);
04520 axis->ResetBit(TAxis::kLabelsUp);
04521 }
04522 Int_t sort = -1;
04523 if (opt.Contains("a")) sort = 0;
04524 if (opt.Contains(">")) sort = 1;
04525 if (opt.Contains("<")) sort = 2;
04526 if (sort < 0) return;
04527 if (sort > 0 && GetDimension() > 2) {
04528 Error("LabelsOption","Sorting by value not implemented for 3-D histograms");
04529 return;
04530 }
04531
04532 Double_t entries = fEntries;
04533 Int_t n = TMath::Min(axis->GetNbins(), labels->GetSize());
04534 Int_t *a = new Int_t[n+2];
04535 Int_t i,j,k;
04536 Double_t *cont = 0;
04537 Double_t *errors = 0;
04538 THashList *labold = new THashList(labels->GetSize(),1);
04539 TIter nextold(labels);
04540 TObject *obj;
04541 while ((obj=nextold())) {
04542 labold->Add(obj);
04543 }
04544 labels->Clear();
04545 if (sort > 0) {
04546
04547 if (GetDimension() == 1) {
04548 cont = new Double_t[n];
04549 if (fSumw2.fN) errors = new Double_t[n];
04550 for (i=1;i<=n;i++) {
04551 cont[i-1] = GetBinContent(i);
04552 if (errors) errors[i-1] = GetBinError(i);
04553 }
04554 if (sort ==1) TMath::Sort(n,cont,a,kTRUE);
04555 else TMath::Sort(n,cont,a,kFALSE);
04556 for (i=1;i<=n;i++) {
04557 SetBinContent(i,cont[a[i-1]]);
04558 if (errors) SetBinError(i,errors[a[i-1]]);
04559 }
04560 for (i=1;i<=n;i++) {
04561 obj = labold->At(a[i-1]);
04562 labels->Add(obj);
04563 obj->SetUniqueID(i);
04564 }
04565 } else if (GetDimension()== 2) {
04566 Double_t *pcont = new Double_t[n+2];
04567 for (i=0;i<=n;i++) pcont[i] = 0;
04568 Int_t nx = fXaxis.GetNbins();
04569 Int_t ny = fYaxis.GetNbins();
04570 cont = new Double_t[(nx+2)*(ny+2)];
04571 if (fSumw2.fN) errors = new Double_t[n];
04572 for (i=1;i<=nx;i++) {
04573 for (j=1;j<=ny;j++) {
04574 cont[i+nx*j] = GetBinContent(i,j);
04575 if (errors) errors[i+nx*j] = GetBinError(i,j);
04576 if (axis == GetXaxis()) k = i;
04577 else k = j;
04578 pcont[k-1] += cont[i+nx*j];
04579 }
04580 }
04581 if (sort ==1) TMath::Sort(n,pcont,a,kTRUE);
04582 else TMath::Sort(n,pcont,a,kFALSE);
04583 for (i=0;i<n;i++) {
04584 obj = labold->At(a[i]);
04585 labels->Add(obj);
04586 obj->SetUniqueID(i+1);
04587 }
04588 delete [] pcont;
04589 for (i=1;i<=nx;i++) {
04590 for (j=1;j<=ny;j++) {
04591 if (axis == GetXaxis()) {
04592 SetBinContent(i,j,cont[a[i-1]+1+nx*j]);
04593 if (errors) SetBinError(i,j,errors[a[i-1]+1+nx*j]);
04594 } else {
04595 SetBinContent(i,j,cont[i+nx*(a[j-1]+1)]);
04596 if (errors) SetBinError(i,j,errors[i+nx*(a[j-1]+1)]);
04597 }
04598 }
04599 }
04600 } else {
04601
04602 }
04603 } else {
04604
04605 const UInt_t kUsed = 1<<18;
04606 TObject *objk=0;
04607 a[0] = 0;
04608 a[n+1] = n+1;
04609 for (i=1;i<=n;i++) {
04610 const char *label = "zzzzzzzzzzzz";
04611 for (j=1;j<=n;j++) {
04612 obj = labold->At(j-1);
04613 if (!obj) continue;
04614 if (obj->TestBit(kUsed)) continue;
04615
04616 if (strcmp(label,obj->GetName()) < 0) continue;
04617 objk = obj;
04618 a[i] = j;
04619 label = obj->GetName();
04620 }
04621 if (objk) {
04622 objk->SetUniqueID(i);
04623 labels->Add(objk);
04624 objk->SetBit(kUsed);
04625 }
04626 }
04627 for (i=1;i<=n;i++) {
04628 obj = labels->At(i-1);
04629 if (!obj) continue;
04630 obj->ResetBit(kUsed);
04631 }
04632
04633 if (GetDimension() == 1) {
04634 cont = new Double_t[n+2];
04635 if (fSumw2.fN) errors = new Double_t[n+2];
04636 for (i=1;i<=n;i++) {
04637 cont[i] = GetBinContent(a[i]);
04638 if (errors) errors[i] = GetBinError(a[i]);
04639 }
04640 for (i=1;i<=n;i++) {
04641 SetBinContent(i,cont[i]);
04642 if (errors) SetBinError(i,errors[i]);
04643 }
04644 } else if (GetDimension()== 2) {
04645 Int_t nx = fXaxis.GetNbins()+2;
04646 Int_t ny = fYaxis.GetNbins()+2;
04647 cont = new Double_t[nx*ny];
04648 if (fSumw2.fN) errors = new Double_t[nx*ny];
04649 for (i=0;i<nx;i++) {
04650 for (j=0;j<ny;j++) {
04651 cont[i+nx*j] = GetBinContent(i,j);
04652 if (errors) errors[i+nx*j] = GetBinError(i,j);
04653 }
04654 }
04655 for (i=0;i<nx;i++) {
04656 for (j=0;j<ny;j++) {
04657 if (axis == GetXaxis()) {
04658 SetBinContent(i,j,cont[a[i]+nx*j]);
04659 if (errors) SetBinError(i,j,errors[a[i]+nx*j]);
04660 } else {
04661 SetBinContent(i,j,cont[i+nx*a[j]]);
04662 if (errors) SetBinError(i,j,errors[i+nx*a[j]]);
04663 }
04664 }
04665 }
04666 } else {
04667 Int_t nx = fXaxis.GetNbins()+2;
04668 Int_t ny = fYaxis.GetNbins()+2;
04669 Int_t nz = fZaxis.GetNbins()+2;
04670 cont = new Double_t[nx*ny*nz];
04671 if (fSumw2.fN) errors = new Double_t[nx*ny*nz];
04672 for (i=0;i<nx;i++) {
04673 for (j=0;j<ny;j++) {
04674 for (k=0;k<nz;k++) {
04675 cont[i+nx*(j+ny*k)] = GetBinContent(i,j,k);
04676 if (errors) errors[i+nx*(j+ny*k)] = GetBinError(i,j,k);
04677 }
04678 }
04679 }
04680 for (i=0;i<nx;i++) {
04681 for (j=0;j<ny;j++) {
04682 for (k=0;k<nz;k++) {
04683 if (axis == GetXaxis()) {
04684 SetBinContent(i,j,k,cont[a[i]+nx*(j+ny*k)]);
04685 if (errors) SetBinError(i,j,k,errors[a[i]+nx*(j+ny*k)]);
04686 } else if (axis == GetYaxis()) {
04687 SetBinContent(i,j,k,cont[i+nx*(a[j]+ny*k)]);
04688 if (errors) SetBinError(i,j,k,errors[i+nx*(a[j]+ny*k)]);
04689 } else {
04690 SetBinContent(i,j,k,cont[i+nx*(j+ny*a[k])]);
04691 if (errors) SetBinError(i,j,k,errors[i+nx*(j+ny*a[k])]);
04692 }
04693 }
04694 }
04695 }
04696 }
04697 }
04698 fEntries = entries;
04699 delete labold;
04700 if (a) delete [] a;
04701 if (cont) delete [] cont;
04702 if (errors) delete [] errors;
04703 }
04704
04705
04706 static Bool_t AlmostEqual(Double_t a, Double_t b, Double_t epsilon = 0.00000001)
04707 {
04708 return TMath::Abs(a - b) < epsilon;
04709 }
04710
04711
04712 static Bool_t AlmostInteger(Double_t a, Double_t epsilon = 0.00000001)
04713 {
04714 return AlmostEqual(a - TMath::Floor(a), 0, epsilon) ||
04715 AlmostEqual(a - TMath::Floor(a), 1, epsilon);
04716 }
04717
04718
04719 Bool_t TH1::SameLimitsAndNBins(const TAxis& axis1, const TAxis& axis2)
04720 {
04721
04722
04723 if ((axis1.GetNbins() == axis2.GetNbins())
04724 && (axis1.GetXmin() == axis2.GetXmin())
04725 && (axis1.GetXmax() == axis2.GetXmax()))
04726 return kTRUE;
04727 else
04728 return kFALSE;
04729 }
04730
04731
04732 Bool_t TH1::RecomputeAxisLimits(TAxis& destAxis, const TAxis& anAxis)
04733 {
04734
04735
04736 if (SameLimitsAndNBins(destAxis, anAxis))
04737 return kTRUE;
04738
04739 if (destAxis.GetXbins()->fN || anAxis.GetXbins()->fN)
04740 return kFALSE;
04741
04742 Double_t width1 = destAxis.GetBinWidth(0);
04743 Double_t width2 = anAxis.GetBinWidth(0);
04744 if (width1 == 0 || width2 == 0)
04745 return kFALSE;
04746
04747 Double_t xmin = TMath::Min(destAxis.GetXmin(), anAxis.GetXmin());
04748 Double_t xmax = TMath::Max(destAxis.GetXmax(), anAxis.GetXmax());
04749 Double_t width = TMath::Max(width1, width2);
04750
04751
04752 if (!AlmostInteger(width/width1) || !AlmostInteger(width/width2))
04753 return kFALSE;
04754
04755
04756 Double_t delta;
04757 delta = (destAxis.GetXmin() - xmin)/width1;
04758 if (!AlmostInteger(delta))
04759 xmin -= (TMath::Ceil(delta) - delta)*width1;
04760
04761 delta = (anAxis.GetXmin() - xmin)/width2;
04762 if (!AlmostInteger(delta))
04763 xmin -= (TMath::Ceil(delta) - delta)*width2;
04764
04765 delta = (destAxis.GetXmin() - xmin)/width1;
04766 if (!AlmostInteger(delta))
04767 return kFALSE;
04768
04769 delta = (xmax - destAxis.GetXmax())/width1;
04770 if (!AlmostInteger(delta))
04771 xmax += (TMath::Ceil(delta) - delta)*width1;
04772
04773 delta = (xmax - anAxis.GetXmax())/width2;
04774 if (!AlmostInteger(delta))
04775 xmax += (TMath::Ceil(delta) - delta)*width2;
04776
04777 delta = (xmax - destAxis.GetXmax())/width1;
04778 if (!AlmostInteger(delta))
04779 return kFALSE;
04780 #ifdef DEBUG
04781 if (!AlmostInteger((xmax - xmin) / width)) {
04782 printf("TH1::RecomputeAxisLimits - Impossible\n");
04783 return kFALSE;
04784 }
04785 #endif
04786 destAxis.Set(TMath::Nint((xmax - xmin)/width), xmin, xmax);
04787 return kTRUE;
04788 }
04789
04790
04791 Long64_t TH1::Merge(TCollection *li)
04792 {
04793
04794
04795
04796
04797
04798
04799
04800
04801
04802
04803
04804
04805
04806
04807
04808
04809
04810
04811
04812
04813
04814
04815
04816
04817
04818
04819
04820
04821
04822
04823
04824
04825
04826
04827
04828
04829 if (!li) return 0;
04830 if (li->IsEmpty()) return (Int_t) GetEntries();
04831
04832
04833
04834 Bool_t mustCleanup = TestBit(kMustCleanup);
04835 if (mustCleanup) ResetBit(kMustCleanup);
04836 TList inlist;
04837 TH1* hclone = (TH1*)Clone("FirstClone");
04838 if (mustCleanup) SetBit(kMustCleanup);
04839 R__ASSERT(hclone);
04840 BufferEmpty(1);
04841 Reset();
04842 SetEntries(0);
04843 inlist.Add(hclone);
04844 inlist.AddAll(li);
04845
04846 THashList allLabels;
04847 THashList* labels=GetXaxis()->GetLabels();
04848 Bool_t haveOneLabel=kFALSE;
04849 if (labels) {
04850 TIter iL(labels);
04851 TObjString* lb;
04852 while ((lb=(TObjString*)iL())) {
04853 haveOneLabel |= (lb && lb->String().Length());
04854 if (!allLabels.FindObject(lb))
04855 allLabels.Add(lb);
04856 }
04857 }
04858
04859 TAxis newXAxis;
04860 Bool_t initialLimitsFound = kFALSE;
04861 Bool_t allHaveLabels = haveOneLabel;
04862 Bool_t same = kTRUE;
04863 Bool_t allHaveLimits = kTRUE;
04864
04865 TIter next(&inlist);
04866 while (TObject *o = next()) {
04867 TH1* h = dynamic_cast<TH1*> (o);
04868 if (!h) {
04869 Error("Add","Attempt to add object of class: %s to a %s",
04870 o->ClassName(),this->ClassName());
04871 return -1;
04872 }
04873 Bool_t hasLimits = h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax();
04874 allHaveLimits = allHaveLimits && hasLimits;
04875
04876 if (hasLimits) {
04877 h->BufferEmpty();
04878 if (!initialLimitsFound) {
04879 initialLimitsFound = kTRUE;
04880 newXAxis.Set(h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
04881 h->GetXaxis()->GetXmax());
04882 }
04883 else {
04884 if (!RecomputeAxisLimits(newXAxis, *(h->GetXaxis()))) {
04885 Error("Merge", "Cannot merge histograms - limits are inconsistent:\n "
04886 "first: (%d, %f, %f), second: (%d, %f, %f)",
04887 newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax(),
04888 h->GetXaxis()->GetNbins(), h->GetXaxis()->GetXmin(),
04889 h->GetXaxis()->GetXmax());
04890 }
04891 }
04892 }
04893 if (allHaveLabels) {
04894 THashList* hlabels=h->GetXaxis()->GetLabels();
04895 Bool_t hasOneLabel=kFALSE;
04896 if (hlabels) {
04897 TIter iL(hlabels);
04898 TObjString* lb;
04899 while ((lb=(TObjString*)iL())) {
04900 hasOneLabel |= (lb && lb->String().Length());
04901 if (!allLabels.FindObject(lb)) {
04902 allLabels.Add(lb);
04903 same = kFALSE;
04904 }
04905 }
04906 }
04907 allHaveLabels&=(labels && hasOneLabel);
04908 if (!allHaveLabels)
04909 Warning("Merge","Not all histograms have labels. I will ignore labels,"
04910 " falling back to bin numbering mode.");
04911 }
04912 }
04913 next.Reset();
04914
04915 same = same && SameLimitsAndNBins(newXAxis, *GetXaxis());
04916 if (!same && initialLimitsFound)
04917 SetBins(newXAxis.GetNbins(), newXAxis.GetXmin(), newXAxis.GetXmax());
04918
04919 if (!allHaveLimits) {
04920
04921 while (TH1* h = (TH1*)next()) {
04922 if (h->GetXaxis()->GetXmin() >= h->GetXaxis()->GetXmax() && h->fBuffer) {
04923
04924 Int_t nbentries = (Int_t)h->fBuffer[0];
04925 for (Int_t i = 0; i < nbentries; i++)
04926 Fill(h->fBuffer[2*i + 2], h->fBuffer[2*i + 1]);
04927
04928
04929 }
04930 }
04931 if (!initialLimitsFound)
04932 return (Int_t) GetEntries();
04933 next.Reset();
04934 }
04935
04936
04937 Double_t stats[kNstat], totstats[kNstat];
04938 for (Int_t i=0;i<kNstat;i++) {totstats[i] = stats[i] = 0;}
04939 GetStats(totstats);
04940 Double_t nentries = GetEntries();
04941 Int_t binx, ix, nx;
04942 Double_t cu;
04943 Bool_t canRebin=TestBit(kCanRebin);
04944 if (!allHaveLabels) ResetBit(kCanRebin);
04945 while (TH1* h=(TH1*)next()) {
04946
04947 if (h->GetXaxis()->GetXmin() < h->GetXaxis()->GetXmax()) {
04948
04949 h->GetStats(stats);
04950 for (Int_t i=0;i<kNstat;i++)
04951 totstats[i] += stats[i];
04952 nentries += h->GetEntries();
04953
04954 nx = h->GetXaxis()->GetNbins();
04955 for (binx = 0; binx <= nx + 1; binx++) {
04956 cu = h->GetBinContent(binx);
04957 if (!allHaveLabels || binx==0 || binx== nx+1) {
04958 if ((!same) && (binx==0 || binx== nx+1)) {
04959 if (cu != 0) {
04960 Error("Merge", "Cannot merge histograms - the histograms have"
04961 " different limits and undeflows/overflows are present."
04962 " The initial histogram is now broken!");
04963 return -1;
04964 }
04965 }
04966 ix = fXaxis.FindBin(h->GetXaxis()->GetBinCenter(binx));
04967 } else {
04968 const char* label=h->GetXaxis()->GetBinLabel(binx);
04969 if (!label || label[0]==0) {
04970 ix = fXaxis.FindBin(h->GetXaxis()->GetBinCenter(binx));
04971 } else {
04972 ix = fXaxis.FindBin(label);
04973 if (ix==-1) ix = fXaxis.FindBin(h->GetXaxis()->GetBinCenter(binx));
04974 }
04975 }
04976 if (ix >= 0) {
04977 AddBinContent(ix,cu);
04978 if (fSumw2.fN) {
04979 Double_t error1 = h->GetBinError(binx);
04980 fSumw2.fArray[ix] += error1*error1;
04981 }
04982 }
04983 }
04984 }
04985 }
04986 if (canRebin) SetBit(kCanRebin);
04987
04988
04989 PutStats(totstats);
04990 SetEntries(nentries);
04991 inlist.Remove(hclone);
04992 delete hclone;
04993 return (Long64_t)nentries;
04994 }
04995
04996
04997 void TH1::Multiply(TF1 *f1, Double_t c1)
04998 {
04999
05000
05001
05002
05003
05004
05005
05006
05007 if (!f1) {
05008 Error("Add","Attempt to multiply by a non-existing function");
05009 return;
05010 }
05011
05012 Int_t nbinsx = GetNbinsX();
05013 Int_t nbinsy = GetNbinsY();
05014 Int_t nbinsz = GetNbinsZ();
05015 if (fDimension < 2) nbinsy = -1;
05016 if (fDimension < 3) nbinsz = -1;
05017
05018
05019 SetMinimum();
05020 SetMaximum();
05021
05022
05023
05024 ResetBit(kCanRebin);
05025
05026
05027 Int_t bin, binx, biny, binz;
05028 Double_t cu,w;
05029 Double_t xx[3];
05030 Double_t *params = 0;
05031 f1->InitArgs(xx,params);
05032 for (binz=0;binz<=nbinsz+1;binz++) {
05033 xx[2] = fZaxis.GetBinCenter(binz);
05034 for (biny=0;biny<=nbinsy+1;biny++) {
05035 xx[1] = fYaxis.GetBinCenter(biny);
05036 for (binx=0;binx<=nbinsx+1;binx++) {
05037 xx[0] = fXaxis.GetBinCenter(binx);
05038 if (!f1->IsInside(xx)) continue;
05039 TF1::RejectPoint(kFALSE);
05040 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
05041 Double_t error1 = GetBinError(bin);
05042 cu = c1*f1->EvalPar(xx);
05043 if (TF1::RejectedPoint()) continue;
05044 w = GetBinContent(bin)*cu;
05045 SetBinContent(bin,w);
05046 if (fSumw2.fN) {
05047 fSumw2.fArray[bin] = cu*cu*error1*error1;
05048 }
05049 }
05050 }
05051 }
05052 ResetStats();
05053 }
05054
05055
05056 void TH1::Multiply(const TH1 *h1)
05057 {
05058
05059
05060
05061
05062
05063
05064
05065
05066
05067
05068
05069
05070
05071 if (!h1) {
05072 Error("Multiply","Attempt to multiply by a non-existing histogram");
05073 return;
05074 }
05075
05076 Int_t nbinsx = GetNbinsX();
05077 Int_t nbinsy = GetNbinsY();
05078 Int_t nbinsz = GetNbinsZ();
05079
05080 try {
05081 CheckConsistency(this,h1);
05082 } catch(DifferentNumberOfBins&) {
05083 Error("Multiply","Attempt to multiply histograms with different number of bins");
05084 return;
05085 } catch(DifferentAxisLimits&) {
05086 Warning("Multiply","Attempt to multiply histograms with different axis limits");
05087 } catch(DifferentBinLimits&) {
05088 Warning("Multiply","Attempt to multiply histograms with different bin limits");
05089 }
05090
05091 if (fDimension < 2) nbinsy = -1;
05092 if (fDimension < 3) nbinsz = -1;
05093
05094
05095 if (fSumw2.fN == 0 && h1->GetSumw2N() != 0) Sumw2();
05096
05097
05098 SetMinimum();
05099 SetMaximum();
05100
05101
05102
05103 ResetBit(kCanRebin);
05104
05105
05106 Int_t bin, binx, biny, binz;
05107 Double_t c0,c1,w;
05108 for (binz=0;binz<=nbinsz+1;binz++) {
05109 for (biny=0;biny<=nbinsy+1;biny++) {
05110 for (binx=0;binx<=nbinsx+1;binx++) {
05111 bin = GetBin(binx,biny,binz);
05112 c0 = GetBinContent(bin);
05113 c1 = h1->GetBinContent(bin);
05114 w = c0*c1;
05115 SetBinContent(bin,w);
05116 if (fSumw2.fN) {
05117 Double_t e0 = GetBinError(bin);
05118 Double_t e1 = h1->GetBinError(bin);
05119 fSumw2.fArray[bin] = (e0*e0*c1*c1 + e1*e1*c0*c0);
05120 }
05121 }
05122 }
05123 }
05124 ResetStats();
05125 }
05126
05127
05128
05129 void TH1::Multiply(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option)
05130 {
05131
05132
05133
05134
05135
05136
05137
05138
05139
05140
05141
05142
05143
05144 TString opt = option;
05145 opt.ToLower();
05146
05147
05148 if (!h1 || !h2) {
05149 Error("Multiply","Attempt to multiply by a non-existing histogram");
05150 return;
05151 }
05152
05153 Int_t nbinsx = GetNbinsX();
05154 Int_t nbinsy = GetNbinsY();
05155 Int_t nbinsz = GetNbinsZ();
05156
05157 try {
05158 CheckConsistency(h1,h2);
05159 CheckConsistency(this,h1);
05160 } catch(DifferentNumberOfBins&) {
05161 Error("Multiply","Attempt to multiply histograms with different number of bins");
05162 return;
05163 } catch(DifferentAxisLimits&) {
05164 Warning("Multiply","Attempt to multiply histograms with different axis limits");
05165 } catch(DifferentBinLimits&) {
05166 Warning("Multiply","Attempt to multiply histograms with different bin limits");
05167 }
05168
05169 if (fDimension < 2) nbinsy = -1;
05170 if (fDimension < 3) nbinsz = -1;
05171
05172
05173 if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
05174
05175
05176 SetMinimum();
05177 SetMaximum();
05178
05179
05180
05181 ResetBit(kCanRebin);
05182
05183
05184 Int_t bin, binx, biny, binz;
05185 Double_t b1,b2,w,d1,d2;
05186 d1 = c1*c1;
05187 d2 = c2*c2;
05188 for (binz=0;binz<=nbinsz+1;binz++) {
05189 for (biny=0;biny<=nbinsy+1;biny++) {
05190 for (binx=0;binx<=nbinsx+1;binx++) {
05191 bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
05192 b1 = h1->GetBinContent(bin);
05193 b2 = h2->GetBinContent(bin);
05194 w = (c1*b1)*(c2*b2);
05195 SetBinContent(bin,w);
05196 if (fSumw2.fN) {
05197 Double_t e1 = h1->GetBinError(bin);
05198 Double_t e2 = h2->GetBinError(bin);
05199 fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1);
05200 }
05201 }
05202 }
05203 }
05204 ResetStats();
05205 }
05206
05207
05208 void TH1::Paint(Option_t *option)
05209 {
05210
05211
05212
05213
05214
05215
05216 GetPainter(option);
05217
05218 if (fPainter) {
05219 if (strlen(option) > 0) fPainter->Paint(option);
05220 else fPainter->Paint(fOption.Data());
05221 }
05222 }
05223
05224
05225 TH1 *TH1::Rebin(Int_t ngroup, const char*newname, const Double_t *xbins)
05226 {
05227
05228
05229
05230
05231
05232
05233
05234
05235
05236
05237
05238
05239
05240
05241
05242
05243
05244
05245
05246
05247
05248
05249
05250
05251
05252
05253
05254
05255
05256
05257
05258
05259
05260
05261 Int_t nbins = fXaxis.GetNbins();
05262 Double_t xmin = fXaxis.GetXmin();
05263 Double_t xmax = fXaxis.GetXmax();
05264 if ((ngroup <= 0) || (ngroup > nbins)) {
05265 Error("Rebin", "Illegal value of ngroup=%d",ngroup);
05266 return 0;
05267 }
05268
05269 if (fDimension > 1 || InheritsFrom(TProfile::Class())) {
05270 Error("Rebin", "Operation valid on 1-D histograms only");
05271 return 0;
05272 }
05273 if (!newname && xbins) {
05274 Error("Rebin","if xbins is specified, newname must be given");
05275 return 0;
05276 }
05277
05278 Int_t newbins = nbins/ngroup;
05279 if (!xbins) {
05280 Int_t nbg = nbins/ngroup;
05281 if (nbg*ngroup != nbins) {
05282 Warning("Rebin", "ngroup=%d must be an exact divider of nbins=%d",ngroup,nbins);
05283 }
05284 }
05285 else {
05286
05287
05288
05289
05290 newbins = ngroup;
05291 ngroup = nbins;
05292 }
05293
05294
05295 Double_t entries = fEntries;
05296 Double_t *oldBins = new Double_t[nbins+2];
05297 Int_t bin, i;
05298 for (bin=0;bin<nbins+2;bin++) oldBins[bin] = GetBinContent(bin);
05299 Double_t *oldErrors = 0;
05300 if (fSumw2.fN != 0) {
05301 oldErrors = new Double_t[nbins+2];
05302 for (bin=0;bin<nbins+2;bin++) oldErrors[bin] = GetBinError(bin);
05303 }
05304
05305
05306 TH1 *hnew = this;
05307 if ((newname && strlen(newname) > 0) || xbins) {
05308 hnew = (TH1*)Clone(newname);
05309 }
05310
05311
05312 Int_t bitRebin = hnew->TestBit(kCanRebin);
05313 hnew->SetBit(kCanRebin,0);
05314
05315
05316 Double_t stat[kNstat];
05317 GetStats(stat);
05318 bool resetStat = false;
05319
05320 if(!xbins && (newbins*ngroup != nbins)) {
05321 xmax = fXaxis.GetBinUpEdge(newbins*ngroup);
05322 resetStat = true;
05323 }
05324
05325 Int_t nDivisions = fXaxis.GetNdivisions();
05326 Color_t axisColor = fXaxis.GetAxisColor();
05327 Color_t labelColor = fXaxis.GetLabelColor();
05328 Style_t labelFont = fXaxis.GetLabelFont();
05329 Float_t labelOffset = fXaxis.GetLabelOffset();
05330 Float_t labelSize = fXaxis.GetLabelSize();
05331 Float_t tickLength = fXaxis.GetTickLength();
05332 Float_t titleOffset = fXaxis.GetTitleOffset();
05333 Float_t titleSize = fXaxis.GetTitleSize();
05334 Color_t titleColor = fXaxis.GetTitleColor();
05335 Style_t titleFont = fXaxis.GetTitleFont();
05336
05337 if(!xbins && (fXaxis.GetXbins()->GetSize() > 0)){
05338 Double_t *bins = new Double_t[newbins+1];
05339 for(i = 0; i <= newbins; ++i) bins[i] = fXaxis.GetBinLowEdge(1+i*ngroup);
05340 hnew->SetBins(newbins,bins);
05341 delete [] bins;
05342 } else if (xbins) {
05343 hnew->SetBins(newbins,xbins);
05344 } else {
05345 hnew->SetBins(newbins,xmin,xmax);
05346 }
05347
05348
05349 fXaxis.SetNdivisions(nDivisions);
05350 fXaxis.SetAxisColor(axisColor);
05351 fXaxis.SetLabelColor(labelColor);
05352 fXaxis.SetLabelFont(labelFont);
05353 fXaxis.SetLabelOffset(labelOffset);
05354 fXaxis.SetLabelSize(labelSize);
05355 fXaxis.SetTickLength(tickLength);
05356 fXaxis.SetTitleOffset(titleOffset);
05357 fXaxis.SetTitleSize(titleSize);
05358 fXaxis.SetTitleColor(titleColor);
05359 fXaxis.SetTitleFont(titleFont);
05360
05361
05362
05363 Int_t startbin = 1;
05364 const Double_t newxmin = hnew->GetXaxis()->GetBinLowEdge(1);
05365 while( fXaxis.GetBinCenter(startbin) < newxmin && startbin <= nbins ) {
05366 startbin++;
05367 }
05368 Int_t oldbin = startbin;
05369 Double_t binContent, binError;
05370 for (bin = 1;bin<=newbins;bin++) {
05371 binContent = 0;
05372 binError = 0;
05373 Int_t imax = ngroup;
05374 Double_t xbinmax = hnew->GetXaxis()->GetBinUpEdge(bin);
05375 for (i=0;i<ngroup;i++) {
05376 if( (hnew == this && (oldbin+i > nbins)) ||
05377 ( hnew != this && (fXaxis.GetBinCenter(oldbin+i) > xbinmax)) ) {
05378 imax = i;
05379 break;
05380 }
05381 binContent += oldBins[oldbin+i];
05382 if (oldErrors) binError += oldErrors[oldbin+i]*oldErrors[oldbin+i];
05383 }
05384 hnew->SetBinContent(bin,binContent);
05385 if (oldErrors) hnew->SetBinError(bin,TMath::Sqrt(binError));
05386 oldbin += imax;
05387 }
05388
05389 binContent = 0;
05390 binError = 0;
05391 for (i = 0; i < startbin; ++i) {
05392 binContent += oldBins[i];
05393 if (oldErrors) binError += oldErrors[i]*oldErrors[i];
05394 }
05395 hnew->SetBinContent(0,binContent);
05396 if (oldErrors) hnew->SetBinError(0,TMath::Sqrt(binError));
05397
05398 binContent = 0;
05399 binError = 0;
05400 for (i = oldbin; i <= nbins+1; ++i) {
05401 binContent += oldBins[i];
05402 if (oldErrors) binError += oldErrors[i]*oldErrors[i];
05403 }
05404 hnew->SetBinContent(newbins+1,binContent);
05405 if (oldErrors) hnew->SetBinError(newbins+1,TMath::Sqrt(binError));
05406 hnew->SetBit(kCanRebin,bitRebin);
05407
05408
05409 hnew->SetEntries(entries);
05410 if (!resetStat) hnew->PutStats(stat);
05411 delete [] oldBins;
05412 if (oldErrors) delete [] oldErrors;
05413 return hnew;
05414 }
05415
05416
05417 Bool_t TH1::FindNewAxisLimits(const TAxis* axis, const Double_t point, Double_t& newMin, Double_t &newMax)
05418 {
05419
05420
05421
05422
05423
05424
05425
05426
05427
05428
05429 Double_t xmin = axis->GetXmin();
05430 Double_t xmax = axis->GetXmax();
05431 if (xmin >= xmax) return kFALSE;
05432 Double_t range = xmax-xmin;
05433 Double_t binsize = range / axis->GetNbins();
05434
05435
05436 Int_t ntimes = 0;
05437 while (point < xmin) {
05438 if (ntimes++ > 64)
05439 return kFALSE;
05440 xmin = xmin - range;
05441 range *= 2;
05442 binsize *= 2;
05443
05444 if ( xmin / binsize - TMath::Floor(xmin / binsize) >= 0.5) {
05445 xmin += 0.5 * binsize;
05446 xmax += 0.5 * binsize;
05447 }
05448 }
05449 while (point >= xmax) {
05450 if (ntimes++ > 64)
05451 return kFALSE;
05452 xmax = xmax + range;
05453 range *= 2;
05454 binsize *= 2;
05455
05456 if ( xmin / binsize - TMath::Floor(xmin / binsize) >= 0.5) {
05457 xmin -= 0.5 * binsize;
05458 xmax -= 0.5 * binsize;
05459 }
05460 }
05461 newMin = xmin;
05462 newMax = xmax;
05463
05464
05465
05466 return kTRUE;
05467 }
05468
05469
05470 void TH1::RebinAxis(Double_t x, TAxis *axis)
05471 {
05472
05473
05474
05475
05476
05477
05478
05479
05480
05481
05482 if (!TestBit(kCanRebin)) return;
05483 if (TMath::IsNaN(x)) {
05484 ResetBit(kCanRebin);
05485 return;
05486 }
05487
05488 if (axis->GetXmin() >= axis->GetXmax()) return;
05489 if (axis->GetNbins() <= 0) return;
05490
05491 Double_t xmin, xmax;
05492 if (!FindNewAxisLimits(axis, x, xmin, xmax))
05493 return;
05494
05495
05496 TH1 *hold = (TH1*)Clone();
05497 hold->SetDirectory(0);
05498
05499 axis->SetLimits(xmin,xmax);
05500
05501 Int_t nbinsx = fXaxis.GetNbins();
05502 Int_t nbinsy = fYaxis.GetNbins();
05503 Int_t nbinsz = fZaxis.GetNbins();
05504
05505
05506 Double_t err,cu;
05507 Double_t bx,by,bz;
05508 Int_t errors = GetSumw2N();
05509 Int_t ix,iy,iz,ibin,binx,biny,binz,bin;
05510 Reset("ICE");
05511 for (binz=1;binz<=nbinsz;binz++) {
05512 bz = hold->GetZaxis()->GetBinCenter(binz);
05513 iz = fZaxis.FindFixBin(bz);
05514 for (biny=1;biny<=nbinsy;biny++) {
05515 by = hold->GetYaxis()->GetBinCenter(biny);
05516 iy = fYaxis.FindFixBin(by);
05517 for (binx=1;binx<=nbinsx;binx++) {
05518 bx = hold->GetXaxis()->GetBinCenter(binx);
05519 ix = fXaxis.FindFixBin(bx);
05520 bin = hold->GetBin(binx,biny,binz);
05521 ibin= GetBin(ix,iy,iz);
05522 cu = hold->GetBinContent(bin);
05523 AddBinContent(ibin,cu);
05524 if (errors) {
05525 err = hold->GetBinError(bin);
05526 fSumw2.fArray[ibin] += err*err;
05527 }
05528 }
05529 }
05530 }
05531 delete hold;
05532 }
05533
05534
05535 void TH1::RecursiveRemove(TObject *obj)
05536 {
05537
05538
05539 if (fFunctions) {
05540 if (!fFunctions->TestBit(kInvalidObject)) fFunctions->RecursiveRemove(obj);
05541 }
05542 }
05543
05544
05545 void TH1::Scale(Double_t c1, Option_t *option)
05546 {
05547
05548
05549
05550
05551
05552
05553
05554
05555
05556
05557
05558
05559
05560
05561
05562
05563
05564
05565
05566 TString opt = option;
05567 opt.ToLower();
05568 Double_t ent = fEntries;
05569 if (opt.Contains("width")) Add(this,this,c1,-1);
05570 else Add(this,this,c1,0);
05571 fEntries = ent;
05572
05573
05574 Int_t ncontours = GetContour();
05575 if (ncontours == 0) return;
05576 Double_t *levels = fContour.GetArray();
05577 for (Int_t i=0;i<ncontours;i++) {
05578 levels[i] *= c1;
05579 }
05580 }
05581
05582
05583
05584 void TH1::SetDefaultBufferSize(Int_t buffersize)
05585 {
05586
05587
05588
05589
05590
05591 if (buffersize < 0) buffersize = 0;
05592 fgBufferSize = buffersize;
05593 }
05594
05595
05596
05597 void TH1::SetDefaultSumw2(Bool_t sumw2)
05598 {
05599
05600
05601
05602
05603
05604 fgDefaultSumw2 = sumw2;
05605 }
05606
05607
05608 void TH1::SetTitle(const char *title)
05609 {
05610
05611
05612
05613
05614
05615
05616
05617
05618 fTitle = title;
05619 fTitle.ReplaceAll("#;",2,"#semicolon",10);
05620
05621
05622 TString str1 = fTitle, str2;
05623 Int_t isc = str1.Index(";");
05624 Int_t lns = str1.Length();
05625
05626 if (isc >=0 ) {
05627 fTitle = str1(0,isc);
05628 str1 = str1(isc+1, lns);
05629 isc = str1.Index(";");
05630 if (isc >=0 ) {
05631 str2 = str1(0,isc);
05632 str2.ReplaceAll("#semicolon",10,";",1);
05633 fXaxis.SetTitle(str2.Data());
05634 lns = str1.Length();
05635 str1 = str1(isc+1, lns);
05636 isc = str1.Index(";");
05637 if (isc >=0 ) {
05638 str2 = str1(0,isc);
05639 str2.ReplaceAll("#semicolon",10,";",1);
05640 fYaxis.SetTitle(str2.Data());
05641 lns = str1.Length();
05642 str1 = str1(isc+1, lns);
05643 str1.ReplaceAll("#semicolon",10,";",1);
05644 fZaxis.SetTitle(str1.Data());
05645 } else {
05646 str1.ReplaceAll("#semicolon",10,";",1);
05647 fYaxis.SetTitle(str1.Data());
05648 }
05649 } else {
05650 str1.ReplaceAll("#semicolon",10,";",1);
05651 fXaxis.SetTitle(str1.Data());
05652 }
05653 }
05654
05655 fTitle.ReplaceAll("#semicolon",10,";",1);
05656
05657 if (gPad && TestBit(kMustCleanup)) gPad->Modified();
05658 }
05659
05660
05661 void TH1::SmoothArray(Int_t nn, Double_t *xx, Int_t ntimes)
05662 {
05663
05664
05665
05666
05667 Int_t ii, jj, ik, jk, kk, nn2;
05668 Double_t hh[6] = {0,0,0,0,0,0};
05669 Double_t *yy = new Double_t[nn];
05670 Double_t *zz = new Double_t[nn];
05671 Double_t *rr = new Double_t[nn];
05672
05673 for (Int_t pass=0;pass<ntimes;pass++) {
05674
05675
05676 for (ii = 0; ii < nn; ii++) {
05677 yy[ii] = xx[ii];
05678 }
05679
05680
05681 for (kk = 1; kk <= 3; kk++) {
05682 ik = 0;
05683 if (kk == 2) ik = 1;
05684 nn2 = nn - ik - 1;
05685
05686
05687 for (ii = ik + 1; ii < nn2; ii++) {
05688 for (jj = 0; jj < 3; jj++) {
05689 hh[jj] = yy[ii + jj - 1];
05690 }
05691 zz[ii] = TMath::Median(3 + 2*ik, hh);
05692 }
05693
05694 if (kk == 1) {
05695
05696 hh[0] = 3*yy[1] - 2*yy[2];
05697 hh[1] = yy[0];
05698 hh[2] = yy[1];
05699 zz[0] = TMath::Median(3, hh);
05700
05701 hh[0] = yy[nn - 2];
05702 hh[1] = yy[nn - 1];
05703 hh[2] = 3*yy[nn - 2] - 2*yy[nn - 3];
05704 zz[nn - 1] = TMath::Median(3, hh);
05705 }
05706 if (kk == 2) {
05707
05708 zz[0] = yy[0];
05709 for (ii = 0; ii < 3; ii++) {
05710 hh[ii] = yy[ii];
05711 }
05712 zz[1] = TMath::Median(3, hh);
05713
05714 for (ii = 0; ii < 3; ii++) {
05715 hh[ii] = yy[nn - 3 + ii];
05716 }
05717 zz[nn - 2] = TMath::Median(3, hh);
05718 zz[nn - 1] = yy[nn - 1];
05719 }
05720 }
05721
05722
05723 for (ii = 2; ii < (nn - 2); ii++) {
05724 if (zz[ii - 1] != zz[ii]) continue;
05725 if (zz[ii] != zz[ii + 1]) continue;
05726 hh[0] = zz[ii - 2] - zz[ii];
05727 hh[1] = zz[ii + 2] - zz[ii];
05728 if (hh[0] * hh[1] < 0) continue;
05729 jk = 1;
05730 if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
05731 yy[ii] = -0.5*zz[ii - 2*jk] + zz[ii]/0.75 + zz[ii + 2*jk] /6.;
05732 yy[ii + jk] = 0.5*(zz[ii + 2*jk] - zz[ii - 2*jk]) + zz[ii];
05733 }
05734
05735
05736 for (ii = 1; ii < nn - 1; ii++) {
05737 rr[ii] = 0.25*yy[ii - 1] + 0.5*yy[ii] + 0.25*yy[ii + 1];
05738 }
05739 rr[0] = yy[0];
05740 rr[nn - 1] = yy[nn - 1];
05741
05742
05743
05744 for (ii = 0; ii < nn; ii++) {
05745 yy[ii] = xx[ii] - rr[ii];
05746 }
05747
05748
05749 for (kk = 1; kk <= 3; kk++) {
05750 ik = 0;
05751 if (kk == 2) ik = 1;
05752 nn2 = nn - ik - 1;
05753
05754
05755 for (ii = ik + 1; ii < nn2; ii++) {
05756 for (jj = 0; jj < 3; jj++) {
05757 hh[jj] = yy[ii + jj - 1];
05758 }
05759 zz[ii] = TMath::Median(3 + 2*ik, hh);
05760 }
05761
05762 if (kk == 1) {
05763
05764 hh[0] = 3*yy[1] - 2*yy[2];
05765 hh[1] = yy[0];
05766 hh[2] = yy[1];
05767 zz[0] = TMath::Median(3, hh);
05768
05769 hh[0] = yy[nn - 2];
05770 hh[1] = yy[nn - 1];
05771 hh[2] = 3*yy[nn - 2] - 2*yy[nn - 3];
05772 zz[nn - 1] = TMath::Median(3, hh);
05773 }
05774 if (kk == 2) {
05775
05776 zz[0] = yy[0];
05777 for (ii = 0; ii < 3; ii++) {
05778 hh[ii] = yy[ii];
05779 }
05780 zz[1] = TMath::Median(3, hh);
05781
05782 for (ii = 0; ii < 3; ii++) {
05783 hh[ii] = yy[nn - 3 + ii];
05784 }
05785 zz[nn - 2] = TMath::Median(3, hh);
05786 zz[nn - 1] = yy[nn - 1];
05787 }
05788 }
05789
05790
05791 for (ii = 2; ii < (nn - 2); ii++) {
05792 if (zz[ii - 1] != zz[ii]) continue;
05793 if (zz[ii] != zz[ii + 1]) continue;
05794 hh[0] = zz[ii - 2] - zz[ii];
05795 hh[1] = zz[ii + 2] - zz[ii];
05796 if (hh[0] * hh[1] < 0) continue;
05797 jk = 1;
05798 if ( TMath::Abs(hh[1]) > TMath::Abs(hh[0]) ) jk = -1;
05799 yy[ii] = -0.5*zz[ii - 2*jk] + zz[ii]/0.75 + zz[ii + 2*jk]/6.;
05800 yy[ii + jk] = 0.5*(zz[ii + 2*jk] - zz[ii - 2*jk]) + zz[ii];
05801 }
05802
05803
05804 for (ii = 1; ii < (nn - 1); ii++) {
05805 zz[ii] = 0.25*yy[ii - 1] + 0.5*yy[ii] + 0.25*yy[ii + 1];
05806 }
05807 zz[0] = yy[0];
05808 zz[nn - 1] = yy[nn - 1];
05809
05810
05811
05812 for (ii = 0; ii < nn; ii++) {
05813 if (xx[ii] < 0) xx[ii] = rr[ii] + zz[ii];
05814 else xx[ii] = TMath::Abs(rr[ii] + zz[ii]);
05815 }
05816 }
05817 delete [] yy;
05818 delete [] zz;
05819 delete [] rr;
05820 }
05821
05822
05823
05824 void TH1::Smooth(Int_t ntimes, Option_t *option)
05825 {
05826
05827
05828
05829
05830
05831
05832
05833 if (fDimension != 1) {
05834 Error("Smooth","Smooth only supported for 1-d histograms");
05835 return;
05836 }
05837 Int_t nbins = fXaxis.GetNbins();
05838 Int_t firstbin = 1, lastbin = nbins;
05839 TString opt = option;
05840 opt.ToLower();
05841 if (opt.Contains("r")) {
05842 firstbin= fXaxis.GetFirst();
05843 lastbin = fXaxis.GetLast();
05844 }
05845 nbins = lastbin - firstbin + 1;
05846 Double_t *xx = new Double_t[nbins];
05847 Double_t nent = fEntries;
05848 Int_t i;
05849 for (i=0;i<nbins;i++) {
05850 xx[i] = GetBinContent(i+firstbin);
05851 }
05852
05853 TH1::SmoothArray(nbins,xx,ntimes);
05854
05855 for (i=0;i<nbins;i++) {
05856 SetBinContent(i+firstbin,xx[i]);
05857 }
05858 fEntries = nent;
05859 delete [] xx;
05860
05861 if (gPad) gPad->Modified();
05862 }
05863
05864
05865
05866 void TH1::StatOverflows(Bool_t flag)
05867 {
05868
05869
05870
05871
05872 fgStatOverflows = flag;
05873 }
05874
05875
05876 void TH1::Streamer(TBuffer &b)
05877 {
05878
05879
05880 if (b.IsReading()) {
05881 UInt_t R__s, R__c;
05882 Version_t R__v = b.ReadVersion(&R__s, &R__c);
05883 if (fDirectory) fDirectory->Remove(this);
05884 fDirectory = 0;
05885 if (R__v > 2) {
05886 b.ReadClassBuffer(TH1::Class(), this, R__v, R__s, R__c);
05887
05888 ResetBit(kMustCleanup);
05889 fXaxis.SetParent(this);
05890 fYaxis.SetParent(this);
05891 fZaxis.SetParent(this);
05892 TIter next(fFunctions);
05893 TObject *obj;
05894 while ((obj=next())) {
05895 if (obj->InheritsFrom(TF1::Class())) ((TF1*)obj)->SetParent(this);
05896 }
05897 return;
05898 }
05899
05900 TNamed::Streamer(b);
05901 TAttLine::Streamer(b);
05902 TAttFill::Streamer(b);
05903 TAttMarker::Streamer(b);
05904 b >> fNcells;
05905 fXaxis.Streamer(b);
05906 fYaxis.Streamer(b);
05907 fZaxis.Streamer(b);
05908 fXaxis.SetParent(this);
05909 fYaxis.SetParent(this);
05910 fZaxis.SetParent(this);
05911 b >> fBarOffset;
05912 b >> fBarWidth;
05913 b >> fEntries;
05914 b >> fTsumw;
05915 b >> fTsumw2;
05916 b >> fTsumwx;
05917 b >> fTsumwx2;
05918 if (R__v < 2) {
05919 Float_t maximum, minimum, norm;
05920 Float_t *contour=0;
05921 b >> maximum; fMaximum = maximum;
05922 b >> minimum; fMinimum = minimum;
05923 b >> norm; fNormFactor = norm;
05924 Int_t n = b.ReadArray(contour);
05925 fContour.Set(n);
05926 for (Int_t i=0;i<n;i++) fContour.fArray[i] = contour[i];
05927 delete [] contour;
05928 } else {
05929 b >> fMaximum;
05930 b >> fMinimum;
05931 b >> fNormFactor;
05932 fContour.Streamer(b);
05933 }
05934 fSumw2.Streamer(b);
05935 fOption.Streamer(b);
05936 fFunctions->Delete();
05937 fFunctions->Streamer(b);
05938 b.CheckByteCount(R__s, R__c, TH1::IsA());
05939
05940 } else {
05941 b.WriteClassBuffer(TH1::Class(),this);
05942 }
05943 }
05944
05945
05946 void TH1::Print(Option_t *option) const
05947 {
05948
05949
05950
05951
05952
05953
05954
05955
05956
05957 printf( "TH1.Print Name = %s, Entries= %d, Total sum= %g\n",GetName(),Int_t(fEntries),GetSumOfWeights());
05958 TString opt = option;
05959 opt.ToLower();
05960 Int_t all;
05961 if (opt.Contains("all")) all = 0;
05962 else if (opt.Contains("range")) all = 1;
05963 else if (opt.Contains("base")) all = 2;
05964 else return;
05965
05966 Int_t bin, binx, biny, binz;
05967 Int_t firstx=0,lastx=0,firsty=0,lasty=0,firstz=0,lastz=0;
05968 if (all == 0) {
05969 lastx = fXaxis.GetNbins()+1;
05970 if (fDimension > 1) lasty = fYaxis.GetNbins()+1;
05971 if (fDimension > 2) lastz = fZaxis.GetNbins()+1;
05972 } else {
05973 firstx = fXaxis.GetFirst(); lastx = fXaxis.GetLast();
05974 if (fDimension > 1) {firsty = fYaxis.GetFirst(); lasty = fYaxis.GetLast();}
05975 if (fDimension > 2) {firstz = fZaxis.GetFirst(); lastz = fZaxis.GetLast();}
05976 }
05977
05978 if (all== 2) {
05979 printf(" Title = %s\n", GetTitle());
05980 printf(" NbinsX= %d, xmin= %g, xmax=%g", fXaxis.GetNbins(), fXaxis.GetXmin(), fXaxis.GetXmax());
05981 if( fDimension > 1) printf(", NbinsY= %d, ymin= %g, ymax=%g", fYaxis.GetNbins(), fYaxis.GetXmin(), fYaxis.GetXmax());
05982 if( fDimension > 2) printf(", NbinsZ= %d, zmin= %g, zmax=%g", fZaxis.GetNbins(), fZaxis.GetXmin(), fZaxis.GetXmax());
05983 printf("\n");
05984 return;
05985 }
05986
05987 Double_t w,e;
05988 Double_t x,y,z;
05989 if (fDimension == 1) {
05990 for (binx=firstx;binx<=lastx;binx++) {
05991 x = fXaxis.GetBinCenter(binx);
05992 w = GetBinContent(binx);
05993 e = GetBinError(binx);
05994 if(fSumw2.fN) printf(" fSumw[%d]=%g, x=%g, error=%g\n",binx,w,x,e);
05995 else printf(" fSumw[%d]=%g, x=%g\n",binx,w,x);
05996 }
05997 }
05998 if (fDimension == 2) {
05999 for (biny=firsty;biny<=lasty;biny++) {
06000 y = fYaxis.GetBinCenter(biny);
06001 for (binx=firstx;binx<=lastx;binx++) {
06002 bin = GetBin(binx,biny);
06003 x = fXaxis.GetBinCenter(binx);
06004 w = GetBinContent(bin);
06005 e = GetBinError(bin);
06006 if(fSumw2.fN) printf(" fSumw[%d][%d]=%g, x=%g, y=%g, error=%g\n",binx,biny,w,x,y,e);
06007 else printf(" fSumw[%d][%d]=%g, x=%g, y=%g\n",binx,biny,w,x,y);
06008 }
06009 }
06010 }
06011 if (fDimension == 3) {
06012 for (binz=firstz;binz<=lastz;binz++) {
06013 z = fZaxis.GetBinCenter(binz);
06014 for (biny=firsty;biny<=lasty;biny++) {
06015 y = fYaxis.GetBinCenter(biny);
06016 for (binx=firstx;binx<=lastx;binx++) {
06017 bin = GetBin(binx,biny,binz);
06018 x = fXaxis.GetBinCenter(binx);
06019 w = GetBinContent(bin);
06020 e = GetBinError(bin);
06021 if(fSumw2.fN) printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%g, error=%g\n",binx,biny,binz,w,x,y,z,e);
06022 else printf(" fSumw[%d][%d][%d]=%g, x=%g, y=%g, z=%g\n",binx,biny,binz,w,x,y,z);
06023 }
06024 }
06025 }
06026 }
06027 }
06028
06029
06030 void TH1::Rebuild(Option_t *)
06031 {
06032
06033
06034 SetBinsLength();
06035 if (fSumw2.fN) {
06036 fSumw2.Set(fNcells);
06037 }
06038 }
06039
06040
06041 void TH1::Reset(Option_t *option)
06042 {
06043
06044
06045
06046
06047
06048
06049 TString opt = option;
06050 opt.ToUpper();
06051 fSumw2.Reset();
06052 if (fIntegral) {delete [] fIntegral; fIntegral = 0;}
06053
06054 if (opt.Contains("M")) {
06055 SetMinimum();
06056 SetMaximum();
06057 }
06058 if (opt.Contains("ICE")) return;
06059 if (fBuffer) {BufferEmpty(); fBuffer[0] = 0;}
06060 fTsumw = 0;
06061 fTsumw2 = 0;
06062 fTsumwx = 0;
06063 fTsumwx2 = 0;
06064 fEntries = 0;
06065
06066 TObject *stats = fFunctions->FindObject("stats");
06067 fFunctions->Remove(stats);
06068
06069
06070
06071
06072 TObject *obj;
06073 while ((obj = fFunctions->First())) {
06074 while(fFunctions->Remove(obj)) { }
06075 delete obj;
06076 }
06077 if(stats) fFunctions->Add(stats);
06078 fContour.Set(0);
06079 }
06080
06081
06082 void TH1::SavePrimitive(ostream &out, Option_t *option )
06083 {
06084
06085
06086 Bool_t nonEqiX = kFALSE;
06087 Bool_t nonEqiY = kFALSE;
06088 Bool_t nonEqiZ = kFALSE;
06089 Int_t i;
06090 static Int_t nxaxis = 0;
06091 static Int_t nyaxis = 0;
06092 static Int_t nzaxis = 0;
06093 TString sxaxis="xAxis",syaxis="yAxis",szaxis="zAxis";
06094
06095
06096
06097 if (GetXaxis()->GetXbins()->fN && GetXaxis()->GetXbins()->fArray) {
06098 nonEqiX = kTRUE;
06099 nxaxis++;
06100 sxaxis += nxaxis;
06101 out << " Double_t "<<sxaxis<<"[" << GetXaxis()->GetXbins()->fN
06102 << "] = {";
06103 for (i = 0; i < GetXaxis()->GetXbins()->fN; i++) {
06104 if (i != 0) out << ", ";
06105 out << GetXaxis()->GetXbins()->fArray[i];
06106 }
06107 out << "}; " << endl;
06108 }
06109
06110
06111
06112 if (fDimension > 1 && GetYaxis()->GetXbins()->fN &&
06113 GetYaxis()->GetXbins()->fArray) {
06114 nonEqiY = kTRUE;
06115 nyaxis++;
06116 syaxis += nyaxis;
06117 out << " Double_t "<<syaxis<<"[" << GetYaxis()->GetXbins()->fN
06118 << "] = {";
06119 for (i = 0; i < GetYaxis()->GetXbins()->fN; i++) {
06120 if (i != 0) out << ", ";
06121 out << GetYaxis()->GetXbins()->fArray[i];
06122 }
06123 out << "}; " << endl;
06124 }
06125
06126
06127
06128 if (fDimension > 2 && GetZaxis()->GetXbins()->fN &&
06129 GetZaxis()->GetXbins()->fArray) {
06130 nonEqiZ = kTRUE;
06131 nzaxis++;
06132 szaxis += nzaxis;
06133 out << " Double_t "<<szaxis<<"[" << GetZaxis()->GetXbins()->fN
06134 << "] = {";
06135 for (i = 0; i < GetZaxis()->GetXbins()->fN; i++) {
06136 if (i != 0) out << ", ";
06137 out << GetZaxis()->GetXbins()->fArray[i];
06138 }
06139 out << "}; " << endl;
06140 }
06141
06142 char quote = '"';
06143 out <<" "<<endl;
06144 out <<" "<< ClassName() <<" *";
06145
06146
06147
06148 static Int_t hcounter = 0;
06149 TString histName = GetName();
06150 if (!fDirectory && !histName.Contains("Graph")) {
06151 hcounter++;
06152 histName += "__";
06153 histName += hcounter;
06154 }
06155 const char *hname = histName.Data();
06156
06157
06158 out << hname << " = new " << ClassName() << "(" << quote
06159 << hname << quote << "," << quote<< GetTitle() << quote
06160 << "," << GetXaxis()->GetNbins();
06161 if (nonEqiX)
06162 out << ", "<<sxaxis;
06163 else
06164 out << "," << GetXaxis()->GetXmin()
06165 << "," << GetXaxis()->GetXmax();
06166 if (fDimension > 1) {
06167 out << "," << GetYaxis()->GetNbins();
06168 if (nonEqiY)
06169 out << ", "<<syaxis;
06170 else
06171 out << "," << GetYaxis()->GetXmin()
06172 << "," << GetYaxis()->GetXmax();
06173 }
06174 if (fDimension > 2) {
06175 out << "," << GetZaxis()->GetNbins();
06176 if (nonEqiZ)
06177 out << ", "<<szaxis;
06178 else
06179 out << "," << GetZaxis()->GetXmin()
06180 << "," << GetZaxis()->GetXmax();
06181 }
06182 out << ");" << endl;
06183
06184
06185 Int_t bin;
06186 for (bin=0;bin<fNcells;bin++) {
06187 Double_t bc = GetBinContent(bin);
06188 if (bc) {
06189 out<<" "<<hname<<"->SetBinContent("<<bin<<","<<bc<<");"<<endl;
06190 }
06191 }
06192
06193
06194 if (fSumw2.fN) {
06195 for (bin=0;bin<fNcells;bin++) {
06196 Double_t be = GetBinError(bin);
06197 if (be) {
06198 out<<" "<<hname<<"->SetBinError("<<bin<<","<<be<<");"<<endl;
06199 }
06200 }
06201 }
06202
06203 TH1::SavePrimitiveHelp(out, hname, option);
06204 }
06205
06206
06207 void TH1::SavePrimitiveHelp(ostream &out, const char *hname, Option_t *option )
06208 {
06209
06210
06211
06212 char quote = '"';
06213 if (TMath::Abs(GetBarOffset()) > 1e-5) {
06214 out<<" "<<hname<<"->SetBarOffset("<<GetBarOffset()<<");"<<endl;
06215 }
06216 if (TMath::Abs(GetBarWidth()-1) > 1e-5) {
06217 out<<" "<<hname<<"->SetBarWidth("<<GetBarWidth()<<");"<<endl;
06218 }
06219 if (fMinimum != -1111) {
06220 out<<" "<<hname<<"->SetMinimum("<<fMinimum<<");"<<endl;
06221 }
06222 if (fMaximum != -1111) {
06223 out<<" "<<hname<<"->SetMaximum("<<fMaximum<<");"<<endl;
06224 }
06225 if (fNormFactor != 0) {
06226 out<<" "<<hname<<"->SetNormFactor("<<fNormFactor<<");"<<endl;
06227 }
06228 if (fEntries != 0) {
06229 out<<" "<<hname<<"->SetEntries("<<fEntries<<");"<<endl;
06230 }
06231 if (fDirectory == 0) {
06232 out<<" "<<hname<<"->SetDirectory(0);"<<endl;
06233 }
06234 if (TestBit(kNoStats)) {
06235 out<<" "<<hname<<"->SetStats(0);"<<endl;
06236 }
06237 if (fOption.Length() != 0) {
06238 out<<" "<<hname<<"->SetOption("<<quote<<fOption.Data()<<quote<<");"<<endl;
06239 }
06240
06241
06242 Int_t ncontours = GetContour();
06243 if (ncontours > 0) {
06244 out<<" "<<hname<<"->SetContour("<<ncontours<<");"<<endl;
06245 Double_t zlevel;
06246 for (Int_t bin=0;bin<ncontours;bin++) {
06247 if (gPad->GetLogz()) {
06248 zlevel = TMath::Power(10,GetContourLevel(bin));
06249 } else {
06250 zlevel = GetContourLevel(bin);
06251 }
06252 out<<" "<<hname<<"->SetContourLevel("<<bin<<","<<zlevel<<");"<<endl;
06253 }
06254 }
06255
06256
06257 TObjOptLink *lnk = (TObjOptLink*)fFunctions->FirstLink();
06258 TObject *obj;
06259 while (lnk) {
06260 obj = lnk->GetObject();
06261 obj->SavePrimitive(out,"nodraw");
06262 if (obj->InheritsFrom(TF1::Class())) {
06263 out<<" "<<hname<<"->GetListOfFunctions()->Add("<<obj->GetName()<<");"<<endl;
06264 } else if (obj->InheritsFrom("TPaveStats")) {
06265 out<<" "<<hname<<"->GetListOfFunctions()->Add(ptstats);"<<endl;
06266 out<<" ptstats->SetParent("<<hname<<"->GetListOfFunctions());"<<endl;
06267 } else {
06268 out<<" "<<hname<<"->GetListOfFunctions()->Add("<<obj->GetName()<<","<<quote<<lnk->GetOption()<<quote<<");"<<endl;
06269 }
06270 lnk = (TObjOptLink*)lnk->Next();
06271 }
06272
06273
06274 SaveFillAttributes(out,hname,0,1001);
06275 SaveLineAttributes(out,hname,1,1,1);
06276 SaveMarkerAttributes(out,hname,1,1,1);
06277 fXaxis.SaveAttributes(out,hname,"->GetXaxis()");
06278 fYaxis.SaveAttributes(out,hname,"->GetYaxis()");
06279 fZaxis.SaveAttributes(out,hname,"->GetZaxis()");
06280 TString opt = option;
06281 opt.ToLower();
06282 if (!opt.Contains("nodraw")) {
06283 out<<" "<<hname<<"->Draw("
06284 <<quote<<option<<quote<<");"<<endl;
06285 }
06286 }
06287
06288
06289 void TH1::UseCurrentStyle()
06290 {
06291
06292
06293 if (!gStyle) return;
06294 if (gStyle->IsReading()) {
06295 fXaxis.ResetAttAxis("X");
06296 fYaxis.ResetAttAxis("Y");
06297 fZaxis.ResetAttAxis("Z");
06298 SetBarOffset(gStyle->GetBarOffset());
06299 SetBarWidth(gStyle->GetBarWidth());
06300 SetFillColor(gStyle->GetHistFillColor());
06301 SetFillStyle(gStyle->GetHistFillStyle());
06302 SetLineColor(gStyle->GetHistLineColor());
06303 SetLineStyle(gStyle->GetHistLineStyle());
06304 SetLineWidth(gStyle->GetHistLineWidth());
06305 SetMarkerColor(gStyle->GetMarkerColor());
06306 SetMarkerStyle(gStyle->GetMarkerStyle());
06307 SetMarkerSize(gStyle->GetMarkerSize());
06308 Int_t dostat = gStyle->GetOptStat();
06309 if (gStyle->GetOptFit() && !dostat) dostat = 1000000001;
06310 SetStats(dostat);
06311 } else {
06312 gStyle->SetBarOffset(fBarOffset);
06313 gStyle->SetBarWidth(fBarWidth);
06314 gStyle->SetHistFillColor(GetFillColor());
06315 gStyle->SetHistFillStyle(GetFillStyle());
06316 gStyle->SetHistLineColor(GetLineColor());
06317 gStyle->SetHistLineStyle(GetLineStyle());
06318 gStyle->SetHistLineWidth(GetLineWidth());
06319 gStyle->SetMarkerColor(GetMarkerColor());
06320 gStyle->SetMarkerStyle(GetMarkerStyle());
06321 gStyle->SetMarkerSize(GetMarkerSize());
06322 gStyle->SetOptStat(TestBit(kNoStats));
06323 }
06324 TIter next(GetListOfFunctions());
06325 TObject *obj;
06326
06327 while ((obj = next())) {
06328 obj->UseCurrentStyle();
06329 }
06330 }
06331
06332
06333 Double_t TH1::GetMean(Int_t axis) const
06334 {
06335
06336
06337
06338
06339
06340
06341
06342
06343
06344
06345
06346
06347
06348
06349
06350
06351
06352
06353
06354
06355
06356
06357
06358 if (axis<1 || (axis>3 && axis<11) || axis>13) return 0;
06359 Double_t stats[kNstat];
06360 for (Int_t i=4;i<kNstat;i++) stats[i] = 0;
06361 GetStats(stats);
06362 if (stats[0] == 0) return 0;
06363 if (axis<4){
06364 Int_t ax[3] = {2,4,7};
06365 return stats[ax[axis-1]]/stats[0];
06366 } else {
06367
06368 Double_t rms = GetRMS(axis-10);
06369 Double_t neff = GetEffectiveEntries();
06370 return ( neff > 0 ? rms/TMath::Sqrt(neff) : 0. );
06371 }
06372 }
06373
06374
06375 Double_t TH1::GetMeanError(Int_t axis) const
06376 {
06377
06378
06379
06380
06381
06382
06383
06384
06385
06386
06387
06388 return GetMean(axis+10);
06389 }
06390
06391
06392 Double_t TH1::GetRMS(Int_t axis) const
06393 {
06394
06395
06396
06397
06398
06399
06400
06401
06402
06403
06404
06405
06406
06407
06408
06409
06410
06411 if (axis<1 || (axis>3 && axis<11) || axis>13) return 0;
06412
06413 Double_t x, rms2, stats[kNstat];
06414 for (Int_t i=4;i<kNstat;i++) stats[i] = 0;
06415 GetStats(stats);
06416 if (stats[0] == 0) return 0;
06417 Int_t ax[3] = {2,4,7};
06418 Int_t axm = ax[axis%10 - 1];
06419 x = stats[axm]/stats[0];
06420 rms2 = TMath::Abs(stats[axm+1]/stats[0] -x*x);
06421 if (axis<10)
06422 return TMath::Sqrt(rms2);
06423 else {
06424
06425
06426 Double_t neff = GetEffectiveEntries();
06427 return ( neff > 0 ? TMath::Sqrt(rms2/(2*neff) ) : 0. );
06428 }
06429 }
06430
06431
06432 Double_t TH1::GetRMSError(Int_t axis) const
06433 {
06434
06435
06436
06437
06438
06439
06440
06441
06442
06443
06444
06445
06446
06447
06448 return GetRMS(axis+10);
06449 }
06450
06451
06452 Double_t TH1::GetSkewness(Int_t axis) const
06453 {
06454
06455
06456
06457
06458
06459
06460
06461 if (axis > 0 && axis <= 3){
06462
06463 Double_t mean = GetMean(axis);
06464 Double_t rms = GetRMS(axis);
06465 Double_t rms3 = rms*rms*rms;
06466
06467 Int_t firstBinX = fXaxis.GetFirst();
06468 Int_t lastBinX = fXaxis.GetLast();
06469 Int_t firstBinY = fYaxis.GetFirst();
06470 Int_t lastBinY = fYaxis.GetLast();
06471 Int_t firstBinZ = fZaxis.GetFirst();
06472 Int_t lastBinZ = fZaxis.GetLast();
06473
06474 if (fgStatOverflows) {
06475 if ( !fXaxis.TestBit(TAxis::kAxisRange) ) {
06476 if (firstBinX == 1) firstBinX = 0;
06477 if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
06478 }
06479 if ( !fYaxis.TestBit(TAxis::kAxisRange) ) {
06480 if (firstBinY == 1) firstBinY = 0;
06481 if (lastBinY == fYaxis.GetNbins() ) lastBinY += 1;
06482 }
06483 if ( !fZaxis.TestBit(TAxis::kAxisRange) ) {
06484 if (firstBinZ == 1) firstBinZ = 0;
06485 if (lastBinZ == fZaxis.GetNbins() ) lastBinZ += 1;
06486 }
06487 }
06488
06489 Double_t x = 0;
06490 Double_t sum=0;
06491 Double_t np=0;
06492 for (Int_t binx = firstBinX; binx <= lastBinX; binx++) {
06493 for (Int_t biny = firstBinY; biny <= lastBinY; biny++) {
06494 for (Int_t binz = firstBinZ; binz <= lastBinZ; binz++) {
06495 if (axis==1 ) x = fXaxis.GetBinCenter(binx);
06496 else if (axis==2 ) x = fYaxis.GetBinCenter(biny);
06497 else if (axis==3 ) x = fZaxis.GetBinCenter(binz);
06498 Double_t w = GetBinContent(binx,biny,binz);
06499 np+=w;
06500 sum+=w*(x-mean)*(x-mean)*(x-mean);
06501 }
06502 }
06503 }
06504 sum/=np*rms3;
06505 return sum;
06506 }
06507 else if (axis > 10 && axis <= 13) {
06508
06509
06510 Double_t neff = GetEffectiveEntries();
06511 return ( neff > 0 ? TMath::Sqrt(6./neff ) : 0. );
06512 }
06513 else {
06514 Error("GetSkewness", "illegal value of parameter");
06515 return 0;
06516 }
06517 }
06518
06519
06520 Double_t TH1::GetKurtosis(Int_t axis) const
06521 {
06522
06523
06524
06525
06526
06527
06528
06529 if (axis > 0 && axis <= 3){
06530
06531 Double_t mean = GetMean(axis);
06532 Double_t rms = GetRMS(axis);
06533 Double_t rms4 = rms*rms*rms*rms;
06534
06535 Int_t firstBinX = fXaxis.GetFirst();
06536 Int_t lastBinX = fXaxis.GetLast();
06537 Int_t firstBinY = fYaxis.GetFirst();
06538 Int_t lastBinY = fYaxis.GetLast();
06539 Int_t firstBinZ = fZaxis.GetFirst();
06540 Int_t lastBinZ = fZaxis.GetLast();
06541
06542 if (fgStatOverflows) {
06543 if ( !fXaxis.TestBit(TAxis::kAxisRange) ) {
06544 if (firstBinX == 1) firstBinX = 0;
06545 if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
06546 }
06547 if ( !fYaxis.TestBit(TAxis::kAxisRange) ) {
06548 if (firstBinY == 1) firstBinY = 0;
06549 if (lastBinY == fYaxis.GetNbins() ) lastBinY += 1;
06550 }
06551 if ( !fZaxis.TestBit(TAxis::kAxisRange) ) {
06552 if (firstBinZ == 1) firstBinZ = 0;
06553 if (lastBinZ == fZaxis.GetNbins() ) lastBinZ += 1;
06554 }
06555 }
06556
06557 Double_t x = 0;
06558 Double_t sum=0;
06559 Double_t np=0;
06560 for (Int_t binx = firstBinX; binx <= lastBinX; binx++) {
06561 for (Int_t biny = firstBinY; biny <= lastBinY; biny++) {
06562 for (Int_t binz = firstBinZ; binz <= lastBinZ; binz++) {
06563 if (axis==1 ) x = fXaxis.GetBinCenter(binx);
06564 else if (axis==2 ) x = fYaxis.GetBinCenter(biny);
06565 else if (axis==3 ) x = fZaxis.GetBinCenter(binz);
06566 Double_t w = GetBinContent(binx,biny,binz);
06567 np+=w;
06568 sum+=w*(x-mean)*(x-mean)*(x-mean)*(x-mean);
06569 }
06570 }
06571 }
06572 sum/=(np*rms4);
06573 return sum-3;
06574
06575 } else if (axis > 10 && axis <= 13) {
06576
06577
06578 Double_t neff = GetEffectiveEntries();
06579 return ( neff > 0 ? TMath::Sqrt(24./neff ) : 0. );
06580 }
06581 else {
06582 Error("GetKurtosis", "illegal value of parameter");
06583 return 0;
06584 }
06585 }
06586
06587
06588
06589 void TH1::GetStats(Double_t *stats) const
06590 {
06591
06592
06593
06594
06595
06596
06597
06598
06599
06600
06601
06602
06603
06604
06605
06606
06607
06608
06609
06610 if (fBuffer) ((TH1*)this)->BufferEmpty();
06611
06612
06613 Int_t bin, binx;
06614 Double_t w,err;
06615 Double_t x;
06616 if ((fTsumw == 0 && fEntries > 0) || fXaxis.TestBit(TAxis::kAxisRange)) {
06617 for (bin=0;bin<4;bin++) stats[bin] = 0;
06618
06619 Int_t firstBinX = fXaxis.GetFirst();
06620 Int_t lastBinX = fXaxis.GetLast();
06621
06622 if (fgStatOverflows && !fXaxis.TestBit(TAxis::kAxisRange)) {
06623 if (firstBinX == 1) firstBinX = 0;
06624 if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
06625 }
06626 for (binx = firstBinX; binx <= lastBinX; binx++) {
06627 x = fXaxis.GetBinCenter(binx);
06628 w = TMath::Abs(GetBinContent(binx));
06629 err = TMath::Abs(GetBinError(binx));
06630 stats[0] += w;
06631 stats[1] += err*err;
06632 stats[2] += w*x;
06633 stats[3] += w*x*x;
06634 }
06635 } else {
06636 stats[0] = fTsumw;
06637 stats[1] = fTsumw2;
06638 stats[2] = fTsumwx;
06639 stats[3] = fTsumwx2;
06640 }
06641 }
06642
06643
06644 void TH1::PutStats(Double_t *stats)
06645 {
06646
06647
06648 fTsumw = stats[0];
06649 fTsumw2 = stats[1];
06650 fTsumwx = stats[2];
06651 fTsumwx2 = stats[3];
06652 }
06653
06654
06655 void TH1::ResetStats()
06656 {
06657
06658
06659
06660
06661 Double_t stats[kNstat] = {0};
06662 fTsumw = 0;
06663 fEntries = 1;
06664 GetStats(stats);
06665 PutStats(stats);
06666 fEntries = TMath::Abs(fTsumw);
06667
06668 if (fSumw2.fN > 0 && fTsumw > 0 && stats[1] > 0 ) fEntries = stats[0]*stats[0]/ stats[1];
06669 }
06670
06671
06672 Double_t TH1::GetSumOfWeights() const
06673 {
06674
06675
06676 Int_t bin,binx,biny,binz;
06677 Double_t sum =0;
06678 for(binz=1; binz<=fZaxis.GetNbins(); binz++) {
06679 for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
06680 for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
06681 bin = GetBin(binx,biny,binz);
06682 sum += GetBinContent(bin);
06683 }
06684 }
06685 }
06686 return sum;
06687 }
06688
06689
06690
06691 Double_t TH1::Integral(Option_t *option) const
06692 {
06693
06694
06695
06696
06697
06698 return Integral(fXaxis.GetFirst(),fXaxis.GetLast(),option);
06699 }
06700
06701
06702 Double_t TH1::Integral(Int_t binx1, Int_t binx2, Option_t *option) const
06703 {
06704
06705
06706
06707
06708 double err = 0;
06709 return DoIntegral(binx1,binx2,0,-1,0,-1,err,option);
06710 }
06711
06712 Double_t TH1::IntegralAndError(Int_t binx1, Int_t binx2, Double_t & error, Option_t *option) const
06713 {
06714
06715
06716
06717
06718
06719
06720 return DoIntegral(binx1,binx2,0,-1,0,-1,error,option,kTRUE);
06721 }
06722
06723
06724
06725 Double_t TH1::DoIntegral(Int_t binx1, Int_t binx2, Int_t biny1, Int_t biny2, Int_t binz1, Int_t binz2, Double_t & error ,
06726 Option_t *option, Bool_t doError) const
06727 {
06728
06729
06730
06731 Int_t nbinsx = GetNbinsX();
06732 if (binx1 < 0) binx1 = 0;
06733 if (binx2 > nbinsx+1 || binx2 < binx1) binx2 = nbinsx+1;
06734 if (GetDimension() > 1) {
06735 Int_t nbinsy = GetNbinsY();
06736 if (biny1 < 0) biny1 = 0;
06737 if (biny2 > nbinsy+1 || biny2 < biny1) biny2 = nbinsy+1;
06738 } else {
06739 biny1 = 0; biny2 = 0;
06740 }
06741 if (GetDimension() > 2) {
06742 Int_t nbinsz = GetNbinsZ();
06743 if (binz1 < 0) binz1 = 0;
06744 if (binz2 > nbinsz+1 || binz2 < binz1) binz2 = nbinsz+1;
06745 } else {
06746 binz1 = 0; binz2 = 0;
06747 }
06748
06749
06750 TString opt = option;
06751 opt.ToLower();
06752 Bool_t width = kFALSE;
06753 if (opt.Contains("width")) width = kTRUE;
06754
06755
06756 Double_t dx = 1.;
06757 Double_t dy = 1.;
06758 Double_t dz = 1.;
06759 Double_t integral = 0;
06760 Double_t igerr2 = 0;
06761 for (Int_t binx = binx1; binx <= binx2; ++binx) {
06762 if (width) dx = fXaxis.GetBinWidth(binx);
06763 for (Int_t biny = biny1; biny <= biny2; ++biny) {
06764 if (width) dy = fYaxis.GetBinWidth(biny);
06765 for (Int_t binz = binz1; binz <= binz2; ++binz) {
06766 if (width) dz = fZaxis.GetBinWidth(binz);
06767 Int_t bin = GetBin(binx, biny, binz);
06768 if (width) integral += GetBinContent(bin)*dx*dy*dz;
06769 else integral += GetBinContent(bin);
06770 if (doError) {
06771 if (width) igerr2 += GetBinError(bin)*GetBinError(bin)*dx*dx*dy*dy*dz*dz;
06772 else igerr2 += GetBinError(bin)*GetBinError(bin);
06773 }
06774 }
06775 }
06776 }
06777
06778 if (doError) error = TMath::Sqrt(igerr2);
06779 return integral;
06780 }
06781
06782
06783
06784 Double_t TH1::KolmogorovTest(const TH1 *h2, Option_t *option) const
06785 {
06786
06787
06788
06789
06790
06791
06792
06793
06794
06795
06796
06797
06798
06799
06800
06801
06802
06803
06804
06805
06806
06807
06808
06809
06810
06811
06812
06813
06814
06815
06816
06817
06818
06819
06820
06821
06822
06823
06824
06825
06826
06827
06828
06829
06830
06831
06832
06833
06834
06835
06836
06837
06838
06839
06840
06841
06842
06843
06844
06845
06846
06847
06848
06849
06850 TString opt = option;
06851 opt.ToUpper();
06852
06853 Double_t prob = 0;
06854 TH1 *h1 = (TH1*)this;
06855 if (h2 == 0) return 0;
06856 TAxis *axis1 = h1->GetXaxis();
06857 TAxis *axis2 = h2->GetXaxis();
06858 Int_t ncx1 = axis1->GetNbins();
06859 Int_t ncx2 = axis2->GetNbins();
06860
06861
06862 if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
06863 Error("KolmogorovTest","Histograms must be 1-D\n");
06864 return 0;
06865 }
06866
06867
06868 if (ncx1 != ncx2) {
06869 Error("KolmogorovTest","Number of channels is different, %d and %d\n",ncx1,ncx2);
06870 return 0;
06871 }
06872
06873
06874 Double_t difprec = 1e-5;
06875 Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
06876 Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
06877 if (diff1 > difprec || diff2 > difprec) {
06878 Error("KolmogorovTest","histograms with different binning");
06879 return 0;
06880 }
06881
06882 Bool_t afunc1 = kFALSE;
06883 Bool_t afunc2 = kFALSE;
06884 Double_t sum1 = 0, sum2 = 0;
06885 Double_t ew1, ew2, w1 = 0, w2 = 0;
06886 Int_t bin;
06887 Int_t ifirst = 1;
06888 Int_t ilast = ncx1;
06889
06890 if (opt.Contains("U")) ifirst = 0;
06891 if (opt.Contains("O")) ilast = ncx1 +1;
06892 for (bin = ifirst; bin <= ilast; bin++) {
06893 sum1 += h1->GetBinContent(bin);
06894 sum2 += h2->GetBinContent(bin);
06895 ew1 = h1->GetBinError(bin);
06896 ew2 = h2->GetBinError(bin);
06897 w1 += ew1*ew1;
06898 w2 += ew2*ew2;
06899 }
06900 if (sum1 == 0) {
06901 Error("KolmogorovTest","Histogram1 %s integral is zero\n",h1->GetName());
06902 return 0;
06903 }
06904 if (sum2 == 0) {
06905 Error("KolmogorovTest","Histogram2 %s integral is zero\n",h2->GetName());
06906 return 0;
06907 }
06908
06909
06910
06911
06912 Double_t esum1 = 0, esum2 = 0;
06913 if (w1 > 0)
06914 esum1 = sum1 * sum1 / w1;
06915 else
06916 afunc1 = kTRUE;
06917
06918 if (w2 > 0)
06919 esum2 = sum2 * sum2 / w2;
06920 else
06921 afunc2 = kTRUE;
06922
06923 if (afunc2 && afunc1) {
06924 Error("KolmogorovTest","Errors are zero for both histograms\n");
06925 return 0;
06926 }
06927
06928
06929 Double_t s1 = 1/sum1;
06930 Double_t s2 = 1/sum2;
06931
06932
06933 Double_t dfmax =0, rsum1 = 0, rsum2 = 0;
06934
06935 for (bin=ifirst;bin<=ilast;bin++) {
06936 rsum1 += s1*h1->GetBinContent(bin);
06937 rsum2 += s2*h2->GetBinContent(bin);
06938 dfmax = TMath::Max(dfmax,TMath::Abs(rsum1-rsum2));
06939 }
06940
06941
06942 Double_t z, prb1=0, prb2=0, prb3=0;
06943
06944
06945 if (afunc1)
06946 z = dfmax*TMath::Sqrt(esum2);
06947
06948 else if (afunc2)
06949 z = dfmax*TMath::Sqrt(esum1);
06950 else
06951
06952 z = dfmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
06953
06954 prob = TMath::KolmogorovProb(z);
06955
06956
06957 if (opt.Contains("N") && !(afunc1 || afunc2 ) ) {
06958
06959 prb1 = prob;
06960 Double_t d12 = esum1-esum2;
06961 Double_t chi2 = d12*d12/(esum1+esum2);
06962 prb2 = TMath::Prob(chi2,1);
06963
06964 if (prob > 0 && prb2 > 0) prob *= prb2*(1-TMath::Log(prob*prb2));
06965 else prob = 0;
06966 }
06967
06968 const Int_t nEXPT = 1000;
06969 if (opt.Contains("X") && !(afunc1 || afunc2 ) ) {
06970 Double_t dSEXPT;
06971 TH1 *hExpt = (TH1*)(gDirectory ? gDirectory->CloneObject(this,kFALSE) : gROOT->CloneObject(this,kFALSE));
06972
06973 prb3 = 0;
06974 for (Int_t i=0; i < nEXPT; i++) {
06975 hExpt->Reset();
06976 hExpt->FillRandom(h1,(Int_t)esum2);
06977 dSEXPT = KolmogorovTest(hExpt,"M");
06978 if (dSEXPT>dfmax) prb3 += 1.0;
06979 }
06980 prb3 /= (Double_t)nEXPT;
06981 delete hExpt;
06982 }
06983
06984
06985 if (opt.Contains("D")) {
06986 printf(" Kolmo Prob h1 = %s, sum bin content =%g effective entries =%g\n",h1->GetName(),sum1,esum1);
06987 printf(" Kolmo Prob h2 = %s, sum bin content =%g effective entries =%g\n",h2->GetName(),sum2,esum2);
06988 printf(" Kolmo Prob = %g, Max Dist = %g\n",prob,dfmax);
06989 if (opt.Contains("N"))
06990 printf(" Kolmo Prob = %f for shape alone, =%f for normalisation alone\n",prb1,prb2);
06991 if (opt.Contains("X"))
06992 printf(" Kolmo Prob = %f with %d pseudo-experiments\n",prb3,nEXPT);
06993 }
06994
06995 if (TMath::Abs(rsum1-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h1=%s\n",h1->GetName());
06996 if (TMath::Abs(rsum2-1) > 0.002) Warning("KolmogorovTest","Numerical problems with h2=%s\n",h2->GetName());
06997
06998 if(opt.Contains("M")) return dfmax;
06999 else if(opt.Contains("X")) return prb3;
07000 else return prob;
07001 }
07002
07003
07004 void TH1::SetContent(const Double_t *content)
07005 {
07006
07007
07008 Int_t bin;
07009 Double_t bincontent;
07010 for (bin=0; bin<fNcells; bin++) {
07011 bincontent = *(content + bin);
07012 SetBinContent(bin, bincontent);
07013 }
07014 }
07015
07016
07017 Int_t TH1::GetContour(Double_t *levels)
07018 {
07019
07020
07021
07022
07023
07024
07025 Int_t nlevels = fContour.fN;
07026 if (levels) {
07027 if (nlevels == 0) {
07028 nlevels = 20;
07029 SetContour(nlevels);
07030 } else {
07031 if (TestBit(kUserContour) == 0) SetContour(nlevels);
07032 }
07033 for (Int_t level=0; level<nlevels; level++) levels[level] = fContour.fArray[level];
07034 }
07035 return nlevels;
07036 }
07037
07038
07039 Double_t TH1::GetContourLevel(Int_t level) const
07040 {
07041
07042
07043
07044 if (level <0 || level >= fContour.fN) return 0;
07045 Double_t zlevel = fContour.fArray[level];
07046 return zlevel;
07047 }
07048
07049
07050 Double_t TH1::GetContourLevelPad(Int_t level) const
07051 {
07052
07053
07054
07055
07056 if (level <0 || level >= fContour.fN) return 0;
07057 Double_t zlevel = fContour.fArray[level];
07058
07059
07060
07061
07062 if (gPad && gPad->GetLogz() && TestBit(kUserContour)) {
07063 if (zlevel <= 0) return 0;
07064 zlevel = TMath::Log10(zlevel);
07065 }
07066 return zlevel;
07067 }
07068
07069
07070 void TH1::SetBuffer(Int_t buffersize, Option_t * )
07071 {
07072
07073
07074 if (fBuffer) {
07075 BufferEmpty();
07076 delete [] fBuffer;
07077 fBuffer = 0;
07078 }
07079 if (buffersize <= 0) {
07080 fBufferSize = 0;
07081 return;
07082 }
07083 if (buffersize < 100) buffersize = 100;
07084 fBufferSize = 1 + buffersize*(fDimension+1);
07085 fBuffer = new Double_t[fBufferSize];
07086 memset(fBuffer,0,8*fBufferSize);
07087 }
07088
07089
07090 void TH1::SetContour(Int_t nlevels, const Double_t *levels)
07091 {
07092
07093
07094
07095
07096
07097
07098
07099
07100 Int_t level;
07101 ResetBit(kUserContour);
07102 if (nlevels <=0 ) {
07103 fContour.Set(0);
07104 return;
07105 }
07106 fContour.Set(nlevels);
07107
07108
07109 if (levels) {
07110 SetBit(kUserContour);
07111 for (level=0; level<nlevels; level++) fContour.fArray[level] = levels[level];
07112 } else {
07113
07114 Double_t zmin = GetMinimum();
07115 Double_t zmax = GetMaximum();
07116 if ((zmin == zmax) && (zmin != 0)) {
07117 zmax += 0.01*TMath::Abs(zmax);
07118 zmin -= 0.01*TMath::Abs(zmin);
07119 }
07120 Double_t dz = (zmax-zmin)/Double_t(nlevels);
07121 if (gPad && gPad->GetLogz()) {
07122 if (zmax <= 0) return;
07123 if (zmin <= 0) zmin = 0.001*zmax;
07124 zmin = TMath::Log10(zmin);
07125 zmax = TMath::Log10(zmax);
07126 dz = (zmax-zmin)/Double_t(nlevels);
07127 }
07128 for (level=0; level<nlevels; level++) {
07129 fContour.fArray[level] = zmin + dz*Double_t(level);
07130 }
07131 }
07132 }
07133
07134
07135 void TH1::SetContourLevel(Int_t level, Double_t value)
07136 {
07137
07138
07139 if (level <0 || level >= fContour.fN) return;
07140 SetBit(kUserContour);
07141 fContour.fArray[level] = value;
07142 }
07143
07144
07145 Double_t TH1::GetMaximum(Double_t maxval) const
07146 {
07147
07148
07149
07150
07151
07152
07153
07154
07155
07156 if (fMaximum != -1111) return fMaximum;
07157 Int_t bin, binx, biny, binz;
07158 Int_t xfirst = fXaxis.GetFirst();
07159 Int_t xlast = fXaxis.GetLast();
07160 Int_t yfirst = fYaxis.GetFirst();
07161 Int_t ylast = fYaxis.GetLast();
07162 Int_t zfirst = fZaxis.GetFirst();
07163 Int_t zlast = fZaxis.GetLast();
07164 Double_t maximum = -FLT_MAX, value;
07165 for (binz=zfirst;binz<=zlast;binz++) {
07166 for (biny=yfirst;biny<=ylast;biny++) {
07167 for (binx=xfirst;binx<=xlast;binx++) {
07168 bin = GetBin(binx,biny,binz);
07169 value = GetBinContent(bin);
07170 if (value > maximum && value < maxval) maximum = value;
07171 }
07172 }
07173 }
07174 return maximum;
07175 }
07176
07177
07178 Int_t TH1::GetMaximumBin() const
07179 {
07180
07181
07182 Int_t locmax, locmay, locmaz;
07183 return GetMaximumBin(locmax, locmay, locmaz);
07184 }
07185
07186
07187 Int_t TH1::GetMaximumBin(Int_t &locmax, Int_t &locmay, Int_t &locmaz) const
07188 {
07189
07190
07191 Int_t bin, binx, biny, binz;
07192 Int_t locm;
07193 Int_t xfirst = fXaxis.GetFirst();
07194 Int_t xlast = fXaxis.GetLast();
07195 Int_t yfirst = fYaxis.GetFirst();
07196 Int_t ylast = fYaxis.GetLast();
07197 Int_t zfirst = fZaxis.GetFirst();
07198 Int_t zlast = fZaxis.GetLast();
07199 Double_t maximum = -FLT_MAX, value;
07200 locm = locmax = locmay = locmaz = 0;
07201 for (binz=zfirst;binz<=zlast;binz++) {
07202 for (biny=yfirst;biny<=ylast;biny++) {
07203 for (binx=xfirst;binx<=xlast;binx++) {
07204 bin = GetBin(binx,biny,binz);
07205 value = GetBinContent(bin);
07206 if (value > maximum) {
07207 maximum = value;
07208 locm = bin;
07209 locmax = binx;
07210 locmay = biny;
07211 locmaz = binz;
07212 }
07213 }
07214 }
07215 }
07216 return locm;
07217 }
07218
07219
07220 Double_t TH1::GetMinimum(Double_t minval) const
07221 {
07222
07223
07224
07225
07226
07227
07228
07229
07230
07231 if (fMinimum != -1111) return fMinimum;
07232 Int_t bin, binx, biny, binz;
07233 Int_t xfirst = fXaxis.GetFirst();
07234 Int_t xlast = fXaxis.GetLast();
07235 Int_t yfirst = fYaxis.GetFirst();
07236 Int_t ylast = fYaxis.GetLast();
07237 Int_t zfirst = fZaxis.GetFirst();
07238 Int_t zlast = fZaxis.GetLast();
07239 Double_t minimum=FLT_MAX, value;
07240 for (binz=zfirst;binz<=zlast;binz++) {
07241 for (biny=yfirst;biny<=ylast;biny++) {
07242 for (binx=xfirst;binx<=xlast;binx++) {
07243 bin = GetBin(binx,biny,binz);
07244 value = GetBinContent(bin);
07245 if (value < minimum && value > minval) minimum = value;
07246 }
07247 }
07248 }
07249 return minimum;
07250 }
07251
07252
07253 Int_t TH1::GetMinimumBin() const
07254 {
07255
07256
07257 Int_t locmix, locmiy, locmiz;
07258 return GetMinimumBin(locmix, locmiy, locmiz);
07259 }
07260
07261
07262 Int_t TH1::GetMinimumBin(Int_t &locmix, Int_t &locmiy, Int_t &locmiz) const
07263 {
07264
07265
07266 Int_t bin, binx, biny, binz;
07267 Int_t locm;
07268 Int_t xfirst = fXaxis.GetFirst();
07269 Int_t xlast = fXaxis.GetLast();
07270 Int_t yfirst = fYaxis.GetFirst();
07271 Int_t ylast = fYaxis.GetLast();
07272 Int_t zfirst = fZaxis.GetFirst();
07273 Int_t zlast = fZaxis.GetLast();
07274 Double_t minimum = FLT_MAX, value;
07275 locm = locmix = locmiy = locmiz = 0;
07276 for (binz=zfirst;binz<=zlast;binz++) {
07277 for (biny=yfirst;biny<=ylast;biny++) {
07278 for (binx=xfirst;binx<=xlast;binx++) {
07279 bin = GetBin(binx,biny,binz);
07280 value = GetBinContent(bin);
07281 if (value < minimum) {
07282 minimum = value;
07283 locm = bin;
07284 locmix = binx;
07285 locmiy = biny;
07286 locmiz = binz;
07287 }
07288 }
07289 }
07290 }
07291 return locm;
07292 }
07293
07294
07295 void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax)
07296 {
07297
07298
07299
07300
07301
07302
07303
07304
07305 if (GetDimension() != 1) {
07306 Error("SetBins","Operation only valid for 1-d histograms");
07307 return;
07308 }
07309 fXaxis.SetRange(0,0);
07310 fXaxis.Set(nx,xmin,xmax);
07311 fYaxis.Set(1,0,1);
07312 fZaxis.Set(1,0,1);
07313 fNcells = nx+2;
07314 SetBinsLength(fNcells);
07315 if (fSumw2.fN) {
07316 fSumw2.Set(fNcells);
07317 }
07318 }
07319
07320
07321 void TH1::SetBins(Int_t nx, const Double_t *xBins)
07322 {
07323
07324
07325
07326
07327
07328
07329
07330
07331 if (GetDimension() != 1) {
07332 Error("SetBins","Operation only valid for 1-d histograms");
07333 return;
07334 }
07335 fXaxis.SetRange(0,0);
07336 fXaxis.Set(nx,xBins);
07337 fYaxis.Set(1,0,1);
07338 fZaxis.Set(1,0,1);
07339 fNcells = nx+2;
07340 SetBinsLength(fNcells);
07341 if (fSumw2.fN) {
07342 fSumw2.Set(fNcells);
07343 }
07344 }
07345
07346
07347 void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax, Int_t ny, Double_t ymin, Double_t ymax)
07348 {
07349
07350
07351
07352
07353
07354
07355
07356
07357 if (GetDimension() != 2) {
07358 Error("SetBins","Operation only valid for 2-D histograms");
07359 return;
07360 }
07361 fXaxis.SetRange(0,0);
07362 fYaxis.SetRange(0,0);
07363 fXaxis.Set(nx,xmin,xmax);
07364 fYaxis.Set(ny,ymin,ymax);
07365 fZaxis.Set(1,0,1);
07366 fNcells = (nx+2)*(ny+2);
07367 SetBinsLength(fNcells);
07368 if (fSumw2.fN) {
07369 fSumw2.Set(fNcells);
07370 }
07371 }
07372
07373
07374 void TH1::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins)
07375 {
07376
07377
07378
07379
07380
07381
07382
07383
07384
07385 if (GetDimension() != 2) {
07386 Error("SetBins","Operation only valid for 2-D histograms");
07387 return;
07388 }
07389 fXaxis.SetRange(0,0);
07390 fYaxis.SetRange(0,0);
07391 fXaxis.Set(nx,xBins);
07392 fYaxis.Set(ny,yBins);
07393 fZaxis.Set(1,0,1);
07394 fNcells = (nx+2)*(ny+2);
07395 SetBinsLength(fNcells);
07396 if (fSumw2.fN) {
07397 fSumw2.Set(fNcells);
07398 }
07399 }
07400
07401
07402
07403 void TH1::SetBins(Int_t nx, Double_t xmin, Double_t xmax, Int_t ny, Double_t ymin, Double_t ymax, Int_t nz, Double_t zmin, Double_t zmax)
07404 {
07405
07406
07407
07408
07409
07410
07411
07412
07413 if (GetDimension() != 3) {
07414 Error("SetBins","Operation only valid for 3-D histograms");
07415 return;
07416 }
07417 fXaxis.SetRange(0,0);
07418 fYaxis.SetRange(0,0);
07419 fZaxis.SetRange(0,0);
07420 fXaxis.Set(nx,xmin,xmax);
07421 fYaxis.Set(ny,ymin,ymax);
07422 fZaxis.Set(nz,zmin,zmax);
07423 fNcells = (nx+2)*(ny+2)*(nz+2);
07424 SetBinsLength(fNcells);
07425 if (fSumw2.fN) {
07426 fSumw2.Set(fNcells);
07427 }
07428 }
07429
07430
07431 void TH1::SetBins(Int_t nx, const Double_t *xBins, Int_t ny, const Double_t *yBins, Int_t nz, const Double_t *zBins)
07432 {
07433
07434
07435
07436
07437
07438
07439
07440
07441
07442
07443 if (GetDimension() != 3) {
07444 Error("SetBins","Operation only valid for 3-D histograms");
07445 return;
07446 }
07447 fXaxis.SetRange(0,0);
07448 fYaxis.SetRange(0,0);
07449 fZaxis.SetRange(0,0);
07450 fXaxis.Set(nx,xBins);
07451 fYaxis.Set(ny,yBins);
07452 fYaxis.Set(nz,zBins);
07453 fNcells = (nx+2)*(ny+2)*(nz+2);
07454 SetBinsLength(fNcells);
07455 if (fSumw2.fN) {
07456 fSumw2.Set(fNcells);
07457 }
07458 }
07459
07460
07461 void TH1::SetMaximum(Double_t maximum)
07462 {
07463
07464
07465
07466
07467
07468
07469
07470
07471
07472
07473
07474 fMaximum = maximum;
07475 }
07476
07477
07478
07479 void TH1::SetMinimum(Double_t minimum)
07480 {
07481
07482
07483
07484
07485
07486
07487
07488
07489
07490
07491
07492 fMinimum = minimum;
07493 }
07494
07495
07496 void TH1::SetDirectory(TDirectory *dir)
07497 {
07498
07499
07500
07501
07502
07503
07504 if (fDirectory == dir) return;
07505 if (fDirectory) fDirectory->Remove(this);
07506 fDirectory = dir;
07507 if (fDirectory) fDirectory->Append(this);
07508 }
07509
07510
07511
07512 void TH1::SetError(const Double_t *error)
07513 {
07514
07515
07516 Int_t bin;
07517 Double_t binerror;
07518 for (bin=0; bin<fNcells; bin++) {
07519 binerror = error[bin];
07520 SetBinError(bin, binerror);
07521 }
07522 }
07523
07524
07525 void TH1::SetName(const char *name)
07526 {
07527
07528
07529
07530
07531
07532 if (fDirectory) fDirectory->Remove(this);
07533 fName = name;
07534 if (fDirectory) fDirectory->Append(this);
07535 }
07536
07537
07538 void TH1::SetNameTitle(const char *name, const char *title)
07539 {
07540
07541
07542
07543
07544
07545 if (fDirectory) fDirectory->Remove(this);
07546 fName = name;
07547 SetTitle(title);
07548 if (fDirectory) fDirectory->Append(this);
07549 }
07550
07551
07552 void TH1::SetStats(Bool_t stats)
07553 {
07554
07555
07556
07557
07558
07559
07560
07561 ResetBit(kNoStats);
07562 if (!stats) {
07563 SetBit(kNoStats);
07564
07565 if (fFunctions) {
07566 TObject *obj = fFunctions->FindObject("stats");
07567 if (obj) {
07568 fFunctions->Remove(obj);
07569 delete obj;
07570 }
07571 }
07572 }
07573 }
07574
07575
07576 void TH1::Sumw2()
07577 {
07578
07579
07580
07581
07582
07583
07584
07585
07586
07587
07588
07589 if (fSumw2.fN == fNcells) {
07590 if (!fgDefaultSumw2 )
07591 Warning("Sumw2","Sum of squares of weights structure already created");
07592 return;
07593 }
07594
07595 fSumw2.Set(fNcells);
07596
07597 if ( fEntries > 0 )
07598 for (Int_t bin=0; bin<fNcells; bin++) {
07599 fSumw2.fArray[bin] = TMath::Abs(GetBinContent(bin));
07600 }
07601 }
07602
07603
07604 TF1 *TH1::GetFunction(const char *name) const
07605 {
07606
07607
07608
07609
07610
07611
07612 return (TF1*)fFunctions->FindObject(name);
07613 }
07614
07615
07616 Double_t TH1::GetBinError(Int_t bin) const
07617 {
07618
07619
07620
07621
07622
07623
07624
07625
07626
07627 if (bin < 0) bin = 0;
07628 if (bin >= fNcells) bin = fNcells-1;
07629 if (fBuffer) ((TH1*)this)->BufferEmpty();
07630 if (fSumw2.fN) {
07631 Double_t err2 = fSumw2.fArray[bin];
07632 return TMath::Sqrt(err2);
07633 }
07634 Double_t error2 = TMath::Abs(GetBinContent(bin));
07635 return TMath::Sqrt(error2);
07636 }
07637
07638
07639 Double_t TH1::GetBinError(Int_t binx, Int_t biny) const
07640 {
07641
07642
07643
07644
07645 Int_t bin = GetBin(binx,biny);
07646 return GetBinError(bin);
07647 }
07648
07649
07650 Double_t TH1::GetBinError(Int_t binx, Int_t biny, Int_t binz) const
07651 {
07652
07653
07654
07655
07656 Int_t bin = GetBin(binx,biny,binz);
07657 return GetBinError(bin);
07658 }
07659
07660
07661 Double_t TH1::GetCellContent(Int_t binx, Int_t biny) const
07662 {
07663
07664
07665
07666
07667 Int_t bin = GetBin(binx,biny);
07668 return GetBinContent(bin);
07669 }
07670
07671
07672 Double_t TH1::GetCellError(Int_t binx, Int_t biny) const
07673 {
07674
07675
07676
07677
07678 Int_t bin = GetBin(binx,biny);
07679 return GetBinError(bin);
07680 }
07681
07682
07683 void TH1::SetBinError(Int_t bin, Double_t error)
07684 {
07685
07686 if (!fSumw2.fN) Sumw2();
07687 if (bin <0 || bin>= fSumw2.fN) return;
07688 fSumw2.fArray[bin] = error*error;
07689 }
07690
07691
07692 void TH1::SetBinContent(Int_t binx, Int_t biny, Double_t content)
07693 {
07694
07695 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07696 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07697 Int_t bin = GetBin(binx,biny);
07698 SetBinContent(bin,content);
07699 }
07700
07701
07702 void TH1::SetBinContent(Int_t binx, Int_t biny, Int_t binz, Double_t content)
07703 {
07704
07705 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07706 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07707 if (binz <0 || binz>fZaxis.GetNbins()+1) return;
07708 Int_t bin = GetBin(binx,biny,binz);
07709 SetBinContent(bin,content);
07710 }
07711
07712
07713 void TH1::SetCellContent(Int_t binx, Int_t biny, Double_t content)
07714 {
07715
07716
07717 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07718 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07719 Int_t bin = GetBin(binx,biny);
07720 SetBinContent(bin,content);
07721 }
07722
07723
07724 void TH1::SetBinError(Int_t binx, Int_t biny, Double_t error)
07725 {
07726
07727 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07728 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07729 Int_t bin = GetBin(binx,biny);
07730 SetBinError(bin,error);
07731 }
07732
07733
07734 void TH1::SetBinError(Int_t binx, Int_t biny, Int_t binz, Double_t error)
07735 {
07736
07737 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07738 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07739 if (binz <0 || binz>fZaxis.GetNbins()+1) return;
07740 Int_t bin = GetBin(binx,biny,binz);
07741 SetBinError(bin,error);
07742 }
07743
07744
07745 void TH1::SetCellError(Int_t binx, Int_t biny, Double_t error)
07746 {
07747
07748 if (binx <0 || binx>fXaxis.GetNbins()+1) return;
07749 if (biny <0 || biny>fYaxis.GetNbins()+1) return;
07750 if (!fSumw2.fN) Sumw2();
07751 Int_t bin = biny*(fXaxis.GetNbins()+2) + binx;
07752 fSumw2.fArray[bin] = error*error;
07753 }
07754
07755
07756 void TH1::SetBinContent(Int_t, Double_t)
07757 {
07758
07759 AbstractMethod("SetBinContent");
07760 }
07761
07762
07763 TH1 *TH1::ShowBackground(Int_t niter, Option_t *option)
07764 {
07765
07766
07767
07768
07769
07770
07771
07772
07773
07774
07775
07776
07777
07778
07779
07780
07781
07782
07783
07784
07785
07786
07787
07788
07789
07790
07791
07792
07793
07794
07795
07796
07797
07798
07799
07800
07801 return (TH1*)gROOT->ProcessLineFast(Form("TSpectrum::StaticBackground((TH1*)0x%lx,%d,\"%s\")",
07802 (ULong_t)this, niter, option));
07803 }
07804
07805
07806 Int_t TH1::ShowPeaks(Double_t sigma, Option_t *option, Double_t threshold)
07807 {
07808
07809
07810
07811
07812
07813
07814
07815 return (Int_t)gROOT->ProcessLineFast(Form("TSpectrum::StaticSearch((TH1*)0x%lx,%g,\"%s\",%g)",
07816 (ULong_t)this, sigma, option, threshold));
07817 }
07818
07819
07820
07821 TH1* TH1::TransformHisto(TVirtualFFT *fft, TH1* h_output, Option_t *option)
07822 {
07823
07824
07825
07826
07827
07828
07829
07830
07831
07832
07833 if (fft->GetNdim()>2){
07834 printf("Only 1d and 2d\n");
07835 return 0;
07836 }
07837 Int_t binx,biny;
07838 TString opt = option;
07839 opt.ToUpper();
07840 Int_t *n = fft->GetN();
07841 TH1 *hout=0;
07842 if (h_output) hout = h_output;
07843 else {
07844 TString name = TString::Format("out_%s", opt.Data());
07845 if (fft->GetNdim()==1)
07846 hout = new TH1D(name, name,n[0], 0, n[0]);
07847 else if (fft->GetNdim()==2)
07848 hout = new TH2D(name, name, n[0], 0, n[0], n[1], 0, n[1]);
07849 }
07850 R__ASSERT(hout != 0);
07851 TString type=fft->GetType();
07852 Int_t ind[2];
07853 if (opt.Contains("RE")){
07854 if (type.Contains("2C") || type.Contains("2HC")) {
07855 Double_t re, im;
07856 for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
07857 for (biny=1; biny<=hout->GetNbinsY(); biny++) {
07858 ind[0] = binx-1; ind[1] = biny-1;
07859 fft->GetPointComplex(ind, re, im);
07860 hout->SetBinContent(binx, biny, re);
07861 }
07862 }
07863 } else {
07864 for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
07865 for (biny=1; biny<=hout->GetNbinsY(); biny++) {
07866 ind[0] = binx-1; ind[1] = biny-1;
07867 hout->SetBinContent(binx, biny, fft->GetPointReal(ind));
07868 }
07869 }
07870 }
07871 }
07872 if (opt.Contains("IM")) {
07873 if (type.Contains("2C") || type.Contains("2HC")) {
07874 Double_t re, im;
07875 for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
07876 for (biny=1; biny<=hout->GetNbinsY(); biny++) {
07877 ind[0] = binx-1; ind[1] = biny-1;
07878 fft->GetPointComplex(ind, re, im);
07879 hout->SetBinContent(binx, biny, im);
07880 }
07881 }
07882 } else {
07883 printf("No complex numbers in the output");
07884 return 0;
07885 }
07886 }
07887 if (opt.Contains("MA")) {
07888 if (type.Contains("2C") || type.Contains("2HC")) {
07889 Double_t re, im;
07890 for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
07891 for (biny=1; biny<=hout->GetNbinsY(); biny++) {
07892 ind[0] = binx-1; ind[1] = biny-1;
07893 fft->GetPointComplex(ind, re, im);
07894 hout->SetBinContent(binx, biny, TMath::Sqrt(re*re + im*im));
07895 }
07896 }
07897 } else {
07898 for (binx = 1; binx<=hout->GetNbinsX(); binx++) {
07899 for (biny=1; biny<=hout->GetNbinsY(); biny++) {
07900 ind[0] = binx-1; ind[1] = biny-1;
07901 hout->SetBinContent(binx, biny, TMath::Abs(fft->GetPointReal(ind)));
07902 }
07903 }
07904 }
07905 }
07906 if (opt.Contains("PH")) {
07907 if (type.Contains("2C") || type.Contains("2HC")){
07908 Double_t re, im, ph;
07909 for (binx = 1; binx<=hout->GetNbinsX(); binx++){
07910 for (biny=1; biny<=hout->GetNbinsY(); biny++){
07911 ind[0] = binx-1; ind[1] = biny-1;
07912 fft->GetPointComplex(ind, re, im);
07913 if (TMath::Abs(re) > 1e-13){
07914 ph = TMath::ATan(im/re);
07915
07916 if (re<0 && im<0)
07917 ph -= TMath::Pi();
07918 if (re<0 && im>=0)
07919 ph += TMath::Pi();
07920 } else {
07921 if (TMath::Abs(im) < 1e-13)
07922 ph = 0;
07923 else if (im>0)
07924 ph = TMath::Pi()*0.5;
07925 else
07926 ph = -TMath::Pi()*0.5;
07927 }
07928 hout->SetBinContent(binx, biny, ph);
07929 }
07930 }
07931 } else {
07932 printf("Pure real output, no phase");
07933 return 0;
07934 }
07935 }
07936
07937 return hout;
07938 }
07939
07940
07941
07942
07943
07944
07945
07946
07947 ClassImp(TH1C)
07948
07949
07950 TH1C::TH1C(): TH1(), TArrayC()
07951 {
07952
07953
07954 fDimension = 1;
07955 SetBinsLength(3);
07956 if (fgDefaultSumw2) Sumw2();
07957 }
07958
07959
07960 TH1C::TH1C(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
07961 : TH1(name,title,nbins,xlow,xup)
07962 {
07963
07964
07965
07966
07967
07968 fDimension = 1;
07969 TArrayC::Set(fNcells);
07970
07971 if (xlow >= xup) SetBuffer(fgBufferSize);
07972 if (fgDefaultSumw2) Sumw2();
07973 }
07974
07975
07976 TH1C::TH1C(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
07977 : TH1(name,title,nbins,xbins)
07978 {
07979
07980
07981
07982
07983
07984 fDimension = 1;
07985 TArrayC::Set(fNcells);
07986 if (fgDefaultSumw2) Sumw2();
07987 }
07988
07989
07990 TH1C::TH1C(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
07991 : TH1(name,title,nbins,xbins)
07992 {
07993
07994
07995
07996
07997
07998 fDimension = 1;
07999 TArrayC::Set(fNcells);
08000 if (fgDefaultSumw2) Sumw2();
08001 }
08002
08003
08004 TH1C::~TH1C()
08005 {
08006
08007 }
08008
08009
08010 TH1C::TH1C(const TH1C &h1c) : TH1(), TArrayC()
08011 {
08012
08013
08014 ((TH1C&)h1c).Copy(*this);
08015 }
08016
08017
08018 void TH1C::AddBinContent(Int_t bin)
08019 {
08020
08021
08022
08023 if (fArray[bin] < 127) fArray[bin]++;
08024 }
08025
08026
08027 void TH1C::AddBinContent(Int_t bin, Double_t w)
08028 {
08029
08030
08031
08032 Int_t newval = fArray[bin] + Int_t(w);
08033 if (newval > -128 && newval < 128) {fArray[bin] = Char_t(newval); return;}
08034 if (newval < -127) fArray[bin] = -127;
08035 if (newval > 127) fArray[bin] = 127;
08036 }
08037
08038
08039 void TH1C::Copy(TObject &newth1) const
08040 {
08041
08042
08043 TH1::Copy(newth1);
08044 }
08045
08046
08047 TH1 *TH1C::DrawCopy(Option_t *option) const
08048 {
08049
08050
08051 TString opt = option;
08052 opt.ToLower();
08053 if (gPad && !opt.Contains("same")) gPad->Clear();
08054 TH1C *newth1 = (TH1C*)Clone();
08055 newth1->SetDirectory(0);
08056 newth1->SetBit(kCanDelete);
08057 newth1->AppendPad(option);
08058 return newth1;
08059 }
08060
08061
08062 Double_t TH1C::GetBinContent(Int_t bin) const
08063 {
08064
08065
08066 if (fBuffer) ((TH1C*)this)->BufferEmpty();
08067 if (bin < 0) bin = 0;
08068 if (bin >= fNcells) bin = fNcells-1;
08069 if (!fArray) return 0;
08070 return Double_t (fArray[bin]);
08071 }
08072
08073
08074 void TH1C::Reset(Option_t *option)
08075 {
08076
08077
08078 TH1::Reset(option);
08079 TArrayC::Reset();
08080 }
08081
08082
08083 void TH1C::SetBinContent(Int_t bin, Double_t content)
08084 {
08085
08086
08087
08088
08089
08090
08091 fEntries++;
08092 fTsumw = 0;
08093 if (bin < 0) return;
08094 if (bin >= fNcells-1) {
08095 if (fXaxis.GetTimeDisplay()) {
08096 while (bin >= fNcells-1) LabelsInflate();
08097 } else {
08098 if (!TestBit(kCanRebin)) {
08099 if (bin == fNcells-1) fArray[bin] = Char_t (content);
08100 return;
08101 }
08102 while (bin >= fNcells-1) LabelsInflate();
08103 }
08104 }
08105 fArray[bin] = Char_t (content);
08106 }
08107
08108
08109 void TH1C::SetBinsLength(Int_t n)
08110 {
08111
08112
08113
08114 if (n < 0) n = fXaxis.GetNbins() + 2;
08115 fNcells = n;
08116 TArrayC::Set(n);
08117 }
08118
08119
08120 TH1C& TH1C::operator=(const TH1C &h1)
08121 {
08122
08123
08124 if (this != &h1) ((TH1C&)h1).Copy(*this);
08125 return *this;
08126 }
08127
08128
08129
08130 TH1C operator*(Double_t c1, const TH1C &h1)
08131 {
08132
08133
08134 TH1C hnew = h1;
08135 hnew.Scale(c1);
08136 hnew.SetDirectory(0);
08137 return hnew;
08138 }
08139
08140
08141 TH1C operator+(const TH1C &h1, const TH1C &h2)
08142 {
08143
08144
08145 TH1C hnew = h1;
08146 hnew.Add(&h2,1);
08147 hnew.SetDirectory(0);
08148 return hnew;
08149 }
08150
08151
08152 TH1C operator-(const TH1C &h1, const TH1C &h2)
08153 {
08154
08155
08156 TH1C hnew = h1;
08157 hnew.Add(&h2,-1);
08158 hnew.SetDirectory(0);
08159 return hnew;
08160 }
08161
08162
08163 TH1C operator*(const TH1C &h1, const TH1C &h2)
08164 {
08165
08166
08167 TH1C hnew = h1;
08168 hnew.Multiply(&h2);
08169 hnew.SetDirectory(0);
08170 return hnew;
08171 }
08172
08173
08174 TH1C operator/(const TH1C &h1, const TH1C &h2)
08175 {
08176
08177
08178 TH1C hnew = h1;
08179 hnew.Divide(&h2);
08180 hnew.SetDirectory(0);
08181 return hnew;
08182 }
08183
08184
08185
08186
08187
08188
08189
08190
08191 ClassImp(TH1S)
08192
08193
08194 TH1S::TH1S(): TH1(), TArrayS()
08195 {
08196
08197
08198 fDimension = 1;
08199 SetBinsLength(3);
08200 if (fgDefaultSumw2) Sumw2();
08201 }
08202
08203
08204 TH1S::TH1S(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
08205 : TH1(name,title,nbins,xlow,xup)
08206 {
08207
08208
08209
08210
08211
08212 fDimension = 1;
08213 TArrayS::Set(fNcells);
08214
08215 if (xlow >= xup) SetBuffer(fgBufferSize);
08216 if (fgDefaultSumw2) Sumw2();
08217 }
08218
08219
08220 TH1S::TH1S(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
08221 : TH1(name,title,nbins,xbins)
08222 {
08223
08224
08225
08226
08227
08228 fDimension = 1;
08229 TArrayS::Set(fNcells);
08230 if (fgDefaultSumw2) Sumw2();
08231 }
08232
08233
08234 TH1S::TH1S(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
08235 : TH1(name,title,nbins,xbins)
08236 {
08237
08238
08239
08240
08241
08242 fDimension = 1;
08243 TArrayS::Set(fNcells);
08244 if (fgDefaultSumw2) Sumw2();
08245 }
08246
08247
08248 TH1S::~TH1S()
08249 {
08250
08251 }
08252
08253
08254 TH1S::TH1S(const TH1S &h1s) : TH1(), TArrayS()
08255 {
08256
08257
08258 ((TH1S&)h1s).Copy(*this);
08259 }
08260
08261
08262 void TH1S::AddBinContent(Int_t bin)
08263 {
08264
08265
08266
08267 if (fArray[bin] < 32767) fArray[bin]++;
08268 }
08269
08270
08271 void TH1S::AddBinContent(Int_t bin, Double_t w)
08272 {
08273
08274
08275
08276 Int_t newval = fArray[bin] + Int_t(w);
08277 if (newval > -32768 && newval < 32768) {fArray[bin] = Short_t(newval); return;}
08278 if (newval < -32767) fArray[bin] = -32767;
08279 if (newval > 32767) fArray[bin] = 32767;
08280 }
08281
08282
08283 void TH1S::Copy(TObject &newth1) const
08284 {
08285
08286
08287 TH1::Copy(newth1);
08288 }
08289
08290
08291 TH1 *TH1S::DrawCopy(Option_t *option) const
08292 {
08293
08294
08295 TString opt = option;
08296 opt.ToLower();
08297 if (gPad && !opt.Contains("same")) gPad->Clear();
08298 TH1S *newth1 = (TH1S*)Clone();
08299 newth1->SetDirectory(0);
08300 newth1->SetBit(kCanDelete);
08301 newth1->AppendPad(option);
08302 return newth1;
08303 }
08304
08305
08306 Double_t TH1S::GetBinContent(Int_t bin) const
08307 {
08308
08309 if (fBuffer) ((TH1S*)this)->BufferEmpty();
08310 if (bin < 0) bin = 0;
08311 if (bin >= fNcells) bin = fNcells-1;
08312 if (!fArray) return 0;
08313 return Double_t (fArray[bin]);
08314 }
08315
08316
08317 void TH1S::Reset(Option_t *option)
08318 {
08319
08320
08321 TH1::Reset(option);
08322 TArrayS::Reset();
08323 }
08324
08325
08326 void TH1S::SetBinContent(Int_t bin, Double_t content)
08327 {
08328
08329
08330
08331
08332
08333
08334 fEntries++;
08335 fTsumw = 0;
08336 if (bin < 0) return;
08337 if (bin >= fNcells-1) {
08338 if (fXaxis.GetTimeDisplay()) {
08339 while (bin >= fNcells-1) LabelsInflate();
08340 } else {
08341 if (!TestBit(kCanRebin)) {
08342 if (bin == fNcells-1) fArray[bin] = Short_t (content);
08343 return;
08344 }
08345 while (bin >= fNcells-1) LabelsInflate();
08346 }
08347 }
08348 fArray[bin] = Short_t (content);
08349 }
08350
08351
08352 void TH1S::SetBinsLength(Int_t n)
08353 {
08354
08355
08356
08357 if (n < 0) n = fXaxis.GetNbins() + 2;
08358 fNcells = n;
08359 TArrayS::Set(n);
08360 }
08361
08362
08363 TH1S& TH1S::operator=(const TH1S &h1)
08364 {
08365
08366
08367 if (this != &h1) ((TH1S&)h1).Copy(*this);
08368 return *this;
08369 }
08370
08371
08372
08373 TH1S operator*(Double_t c1, const TH1S &h1)
08374 {
08375
08376
08377 TH1S hnew = h1;
08378 hnew.Scale(c1);
08379 hnew.SetDirectory(0);
08380 return hnew;
08381 }
08382
08383
08384 TH1S operator+(const TH1S &h1, const TH1S &h2)
08385 {
08386
08387
08388 TH1S hnew = h1;
08389 hnew.Add(&h2,1);
08390 hnew.SetDirectory(0);
08391 return hnew;
08392 }
08393
08394
08395 TH1S operator-(const TH1S &h1, const TH1S &h2)
08396 {
08397
08398
08399 TH1S hnew = h1;
08400 hnew.Add(&h2,-1);
08401 hnew.SetDirectory(0);
08402 return hnew;
08403 }
08404
08405
08406 TH1S operator*(const TH1S &h1, const TH1S &h2)
08407 {
08408
08409
08410 TH1S hnew = h1;
08411 hnew.Multiply(&h2);
08412 hnew.SetDirectory(0);
08413 return hnew;
08414 }
08415
08416
08417 TH1S operator/(const TH1S &h1, const TH1S &h2)
08418 {
08419
08420
08421 TH1S hnew = h1;
08422 hnew.Divide(&h2);
08423 hnew.SetDirectory(0);
08424 return hnew;
08425 }
08426
08427
08428
08429
08430
08431
08432
08433 ClassImp(TH1I)
08434
08435
08436 TH1I::TH1I(): TH1(), TArrayI()
08437 {
08438
08439
08440 fDimension = 1;
08441 SetBinsLength(3);
08442 if (fgDefaultSumw2) Sumw2();
08443 }
08444
08445
08446 TH1I::TH1I(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
08447 : TH1(name,title,nbins,xlow,xup)
08448 {
08449
08450
08451
08452
08453
08454 fDimension = 1;
08455 TArrayI::Set(fNcells);
08456
08457 if (xlow >= xup) SetBuffer(fgBufferSize);
08458 if (fgDefaultSumw2) Sumw2();
08459 }
08460
08461
08462 TH1I::TH1I(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
08463 : TH1(name,title,nbins,xbins)
08464 {
08465
08466
08467
08468
08469
08470 fDimension = 1;
08471 TArrayI::Set(fNcells);
08472 if (fgDefaultSumw2) Sumw2();
08473 }
08474
08475
08476 TH1I::TH1I(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
08477 : TH1(name,title,nbins,xbins)
08478 {
08479
08480
08481
08482
08483
08484 fDimension = 1;
08485 TArrayI::Set(fNcells);
08486 if (fgDefaultSumw2) Sumw2();
08487 }
08488
08489
08490 TH1I::~TH1I()
08491 {
08492
08493 }
08494
08495
08496 TH1I::TH1I(const TH1I &h1i) : TH1(), TArrayI()
08497 {
08498
08499
08500 ((TH1I&)h1i).Copy(*this);
08501 }
08502
08503
08504 void TH1I::AddBinContent(Int_t bin)
08505 {
08506
08507
08508
08509 if (fArray[bin] < 2147483647) fArray[bin]++;
08510 }
08511
08512
08513 void TH1I::AddBinContent(Int_t bin, Double_t w)
08514 {
08515
08516
08517
08518 Int_t newval = fArray[bin] + Int_t(w);
08519 if (newval > -2147483647 && newval < 2147483647) {fArray[bin] = Int_t(newval); return;}
08520 if (newval < -2147483647) fArray[bin] = -2147483647;
08521 if (newval > 2147483647) fArray[bin] = 2147483647;
08522 }
08523
08524
08525 void TH1I::Copy(TObject &newth1) const
08526 {
08527
08528
08529 TH1::Copy(newth1);
08530 }
08531
08532
08533 TH1 *TH1I::DrawCopy(Option_t *option) const
08534 {
08535
08536
08537 TString opt = option;
08538 opt.ToLower();
08539 if (gPad && !opt.Contains("same")) gPad->Clear();
08540 TH1I *newth1 = (TH1I*)Clone();
08541 newth1->SetDirectory(0);
08542 newth1->SetBit(kCanDelete);
08543 newth1->AppendPad(option);
08544 return newth1;
08545 }
08546
08547
08548 Double_t TH1I::GetBinContent(Int_t bin) const
08549 {
08550
08551 if (fBuffer) ((TH1I*)this)->BufferEmpty();
08552 if (bin < 0) bin = 0;
08553 if (bin >= fNcells) bin = fNcells-1;
08554 if (!fArray) return 0;
08555 return Double_t (fArray[bin]);
08556 }
08557
08558
08559 void TH1I::Reset(Option_t *option)
08560 {
08561
08562
08563 TH1::Reset(option);
08564 TArrayI::Reset();
08565 }
08566
08567
08568 void TH1I::SetBinContent(Int_t bin, Double_t content)
08569 {
08570
08571
08572
08573
08574
08575
08576 fEntries++;
08577 fTsumw = 0;
08578 if (bin < 0) return;
08579 if (bin >= fNcells-1) {
08580 if (fXaxis.GetTimeDisplay()) {
08581 while (bin >= fNcells-1) LabelsInflate();
08582 } else {
08583 if (!TestBit(kCanRebin)) {
08584 if (bin == fNcells-1) fArray[bin] = Int_t (content);
08585 return;
08586 }
08587 while (bin >= fNcells-1) LabelsInflate();
08588 }
08589 }
08590 fArray[bin] = Int_t (content);
08591 }
08592
08593
08594 void TH1I::SetBinsLength(Int_t n)
08595 {
08596
08597
08598
08599 if (n < 0) n = fXaxis.GetNbins() + 2;
08600 fNcells = n;
08601 TArrayI::Set(n);
08602 }
08603
08604
08605 TH1I& TH1I::operator=(const TH1I &h1)
08606 {
08607
08608
08609 if (this != &h1) ((TH1I&)h1).Copy(*this);
08610 return *this;
08611 }
08612
08613
08614
08615 TH1I operator*(Double_t c1, const TH1I &h1)
08616 {
08617
08618
08619 TH1I hnew = h1;
08620 hnew.Scale(c1);
08621 hnew.SetDirectory(0);
08622 return hnew;
08623 }
08624
08625
08626 TH1I operator+(const TH1I &h1, const TH1I &h2)
08627 {
08628
08629
08630 TH1I hnew = h1;
08631 hnew.Add(&h2,1);
08632 hnew.SetDirectory(0);
08633 return hnew;
08634 }
08635
08636
08637 TH1I operator-(const TH1I &h1, const TH1I &h2)
08638 {
08639
08640
08641 TH1I hnew = h1;
08642 hnew.Add(&h2,-1);
08643 hnew.SetDirectory(0);
08644 return hnew;
08645 }
08646
08647
08648 TH1I operator*(const TH1I &h1, const TH1I &h2)
08649 {
08650
08651
08652 TH1I hnew = h1;
08653 hnew.Multiply(&h2);
08654 hnew.SetDirectory(0);
08655 return hnew;
08656 }
08657
08658
08659 TH1I operator/(const TH1I &h1, const TH1I &h2)
08660 {
08661
08662
08663 TH1I hnew = h1;
08664 hnew.Divide(&h2);
08665 hnew.SetDirectory(0);
08666 return hnew;
08667 }
08668
08669
08670
08671
08672
08673
08674
08675 ClassImp(TH1F)
08676
08677
08678 TH1F::TH1F(): TH1(), TArrayF()
08679 {
08680
08681
08682 fDimension = 1;
08683 SetBinsLength(3);
08684 if (fgDefaultSumw2) Sumw2();
08685 }
08686
08687
08688 TH1F::TH1F(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
08689 : TH1(name,title,nbins,xlow,xup)
08690 {
08691
08692
08693
08694
08695
08696 fDimension = 1;
08697 TArrayF::Set(fNcells);
08698
08699 if (xlow >= xup) SetBuffer(fgBufferSize);
08700 if (fgDefaultSumw2) Sumw2();
08701 }
08702
08703
08704 TH1F::TH1F(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
08705 : TH1(name,title,nbins,xbins)
08706 {
08707
08708
08709
08710
08711
08712 fDimension = 1;
08713 TArrayF::Set(fNcells);
08714 if (fgDefaultSumw2) Sumw2();
08715 }
08716
08717
08718 TH1F::TH1F(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
08719 : TH1(name,title,nbins,xbins)
08720 {
08721
08722
08723
08724
08725
08726 fDimension = 1;
08727 TArrayF::Set(fNcells);
08728 if (fgDefaultSumw2) Sumw2();
08729 }
08730
08731
08732 TH1F::TH1F(const TVectorF &v)
08733 : TH1("TVectorF","",v.GetNrows(),0,v.GetNrows())
08734 {
08735
08736
08737
08738 TArrayF::Set(fNcells);
08739 fDimension = 1;
08740 Int_t ivlow = v.GetLwb();
08741 for (Int_t i=0;i<fNcells-2;i++) {
08742 SetBinContent(i+1,v(i+ivlow));
08743 }
08744 TArrayF::Set(fNcells);
08745 if (fgDefaultSumw2) Sumw2();
08746 }
08747
08748
08749 TH1F::TH1F(const TH1F &h) : TH1(), TArrayF()
08750 {
08751
08752
08753 ((TH1F&)h).Copy(*this);
08754 }
08755
08756
08757 TH1F::~TH1F()
08758 {
08759
08760 }
08761
08762
08763 void TH1F::Copy(TObject &newth1) const
08764 {
08765
08766
08767 TH1::Copy(newth1);
08768 }
08769
08770
08771 TH1 *TH1F::DrawCopy(Option_t *option) const
08772 {
08773
08774
08775 TString opt = option;
08776 opt.ToLower();
08777 if (gPad && !opt.Contains("same")) gPad->Clear();
08778 TH1F *newth1 = (TH1F*)Clone();
08779 newth1->SetDirectory(0);
08780 newth1->SetBit(kCanDelete);
08781 newth1->AppendPad(option);
08782 return newth1;
08783 }
08784
08785
08786 Double_t TH1F::GetBinContent(Int_t bin) const
08787 {
08788
08789 if (fBuffer) ((TH1F*)this)->BufferEmpty();
08790 if (bin < 0) bin = 0;
08791 if (bin >= fNcells) bin = fNcells-1;
08792 if (!fArray) return 0;
08793 return Double_t (fArray[bin]);
08794 }
08795
08796
08797 void TH1F::Reset(Option_t *option)
08798 {
08799
08800
08801 TH1::Reset(option);
08802 TArrayF::Reset();
08803 }
08804
08805
08806 void TH1F::SetBinContent(Int_t bin, Double_t content)
08807 {
08808
08809
08810
08811
08812
08813
08814 fEntries++;
08815 fTsumw = 0;
08816 if (bin < 0) return;
08817 if (bin >= fNcells-1) {
08818 if (fXaxis.GetTimeDisplay()) {
08819 while (bin >= fNcells-1) LabelsInflate();
08820 } else {
08821 if (!TestBit(kCanRebin)) {
08822 if (bin == fNcells-1) fArray[bin] = Float_t (content);
08823 return;
08824 }
08825 while (bin >= fNcells-1) LabelsInflate();
08826 }
08827 }
08828 fArray[bin] = Float_t (content);
08829 }
08830
08831
08832 void TH1F::SetBinsLength(Int_t n)
08833 {
08834
08835
08836
08837 if (n < 0) n = fXaxis.GetNbins() + 2;
08838 fNcells = n;
08839 TArrayF::Set(n);
08840 }
08841
08842
08843 TH1F& TH1F::operator=(const TH1F &h1)
08844 {
08845
08846
08847 if (this != &h1) ((TH1F&)h1).Copy(*this);
08848 return *this;
08849 }
08850
08851
08852
08853 TH1F operator*(Double_t c1, const TH1F &h1)
08854 {
08855
08856
08857 TH1F hnew = h1;
08858 hnew.Scale(c1);
08859 hnew.SetDirectory(0);
08860 return hnew;
08861 }
08862
08863
08864 TH1F operator+(const TH1F &h1, const TH1F &h2)
08865 {
08866
08867
08868 TH1F hnew = h1;
08869 hnew.Add(&h2,1);
08870 hnew.SetDirectory(0);
08871 return hnew;
08872 }
08873
08874
08875 TH1F operator-(const TH1F &h1, const TH1F &h2)
08876 {
08877
08878
08879 TH1F hnew = h1;
08880 hnew.Add(&h2,-1);
08881 hnew.SetDirectory(0);
08882 return hnew;
08883 }
08884
08885
08886 TH1F operator*(const TH1F &h1, const TH1F &h2)
08887 {
08888
08889
08890 TH1F hnew = h1;
08891 hnew.Multiply(&h2);
08892 hnew.SetDirectory(0);
08893 return hnew;
08894 }
08895
08896
08897 TH1F operator/(const TH1F &h1, const TH1F &h2)
08898 {
08899
08900
08901 TH1F hnew = h1;
08902 hnew.Divide(&h2);
08903 hnew.SetDirectory(0);
08904 return hnew;
08905 }
08906
08907
08908
08909
08910
08911
08912
08913
08914 ClassImp(TH1D)
08915
08916
08917 TH1D::TH1D(): TH1(), TArrayD()
08918 {
08919
08920
08921 fDimension = 1;
08922 SetBinsLength(3);
08923 if (fgDefaultSumw2) Sumw2();
08924 }
08925
08926
08927 TH1D::TH1D(const char *name,const char *title,Int_t nbins,Double_t xlow,Double_t xup)
08928 : TH1(name,title,nbins,xlow,xup)
08929 {
08930
08931
08932
08933
08934
08935 fDimension = 1;
08936 TArrayD::Set(fNcells);
08937
08938 if (xlow >= xup) SetBuffer(fgBufferSize);
08939 if (fgDefaultSumw2) Sumw2();
08940 }
08941
08942
08943 TH1D::TH1D(const char *name,const char *title,Int_t nbins,const Float_t *xbins)
08944 : TH1(name,title,nbins,xbins)
08945 {
08946
08947
08948
08949
08950
08951 fDimension = 1;
08952 TArrayD::Set(fNcells);
08953 if (fgDefaultSumw2) Sumw2();
08954 }
08955
08956
08957 TH1D::TH1D(const char *name,const char *title,Int_t nbins,const Double_t *xbins)
08958 : TH1(name,title,nbins,xbins)
08959 {
08960
08961
08962
08963
08964
08965 fDimension = 1;
08966 TArrayD::Set(fNcells);
08967 if (fgDefaultSumw2) Sumw2();
08968 }
08969
08970
08971 TH1D::TH1D(const TVectorD &v)
08972 : TH1("TVectorD","",v.GetNrows(),0,v.GetNrows())
08973 {
08974
08975
08976
08977 TArrayD::Set(fNcells);
08978 fDimension = 1;
08979 Int_t ivlow = v.GetLwb();
08980 for (Int_t i=0;i<fNcells-2;i++) {
08981 SetBinContent(i+1,v(i+ivlow));
08982 }
08983 TArrayD::Set(fNcells);
08984 if (fgDefaultSumw2) Sumw2();
08985 }
08986
08987
08988 TH1D::~TH1D()
08989 {
08990
08991 }
08992
08993
08994 TH1D::TH1D(const TH1D &h1d) : TH1(), TArrayD()
08995 {
08996
08997
08998 ((TH1D&)h1d).Copy(*this);
08999 }
09000
09001
09002 void TH1D::Copy(TObject &newth1) const
09003 {
09004
09005
09006 TH1::Copy(newth1);
09007 }
09008
09009
09010 TH1 *TH1D::DrawCopy(Option_t *option) const
09011 {
09012
09013
09014 TString opt = option;
09015 opt.ToLower();
09016 if (gPad && !opt.Contains("same")) gPad->Clear();
09017 TH1D *newth1 = (TH1D*)Clone();
09018 newth1->SetDirectory(0);
09019 newth1->SetBit(kCanDelete);
09020 newth1->AppendPad(option);
09021 return newth1;
09022 }
09023
09024
09025 Double_t TH1D::GetBinContent(Int_t bin) const
09026 {
09027
09028 if (fBuffer) ((TH1D*)this)->BufferEmpty();
09029 if (bin < 0) bin = 0;
09030 if (bin >= fNcells) bin = fNcells-1;
09031 if (!fArray) return 0;
09032 return Double_t (fArray[bin]);
09033 }
09034
09035
09036 void TH1D::Reset(Option_t *option)
09037 {
09038
09039
09040 TH1::Reset(option);
09041 TArrayD::Reset();
09042 }
09043
09044
09045 void TH1D::SetBinContent(Int_t bin, Double_t content)
09046 {
09047
09048
09049
09050
09051
09052
09053 fEntries++;
09054 fTsumw = 0;
09055 if (bin < 0) return;
09056 if (bin >= fNcells-1) {
09057 if (fXaxis.GetTimeDisplay()) {
09058 while (bin >= fNcells-1) LabelsInflate();
09059 } else {
09060 if (!TestBit(kCanRebin)) {
09061 if (bin == fNcells-1) fArray[bin] = content;
09062 return;
09063 }
09064 while (bin >= fNcells-1) LabelsInflate();
09065 }
09066 }
09067 fArray[bin] = content;
09068 }
09069
09070
09071 void TH1D::SetBinsLength(Int_t n)
09072 {
09073
09074
09075
09076 if (n < 0) n = fXaxis.GetNbins() + 2;
09077 fNcells = n;
09078 TArrayD::Set(n);
09079 }
09080
09081
09082 TH1D& TH1D::operator=(const TH1D &h1)
09083 {
09084
09085
09086 if (this != &h1) ((TH1D&)h1).Copy(*this);
09087 return *this;
09088 }
09089
09090
09091 TH1D operator*(Double_t c1, const TH1D &h1)
09092 {
09093
09094
09095 TH1D hnew = h1;
09096 hnew.Scale(c1);
09097 hnew.SetDirectory(0);
09098 return hnew;
09099 }
09100
09101
09102 TH1D operator+(const TH1D &h1, const TH1D &h2)
09103 {
09104
09105
09106 TH1D hnew = h1;
09107 hnew.Add(&h2,1);
09108 hnew.SetDirectory(0);
09109 return hnew;
09110 }
09111
09112
09113 TH1D operator-(const TH1D &h1, const TH1D &h2)
09114 {
09115
09116
09117 TH1D hnew = h1;
09118 hnew.Add(&h2,-1);
09119 hnew.SetDirectory(0);
09120 return hnew;
09121 }
09122
09123
09124 TH1D operator*(const TH1D &h1, const TH1D &h2)
09125 {
09126
09127
09128 TH1D hnew = h1;
09129 hnew.Multiply(&h2);
09130 hnew.SetDirectory(0);
09131 return hnew;
09132 }
09133
09134
09135 TH1D operator/(const TH1D &h1, const TH1D &h2)
09136 {
09137
09138
09139 TH1D hnew = h1;
09140 hnew.Divide(&h2);
09141 hnew.SetDirectory(0);
09142 return hnew;
09143 }
09144
09145
09146 TH1 *R__H(Int_t hid)
09147 {
09148
09149
09150
09151
09152 TString hname;
09153 if(hid >= 0) hname.Form("h%d",hid);
09154 else hname.Form("h_%d",hid);
09155 return (TH1*)gDirectory->Get(hname);
09156 }
09157
09158
09159 TH1 *R__H(const char * hname)
09160 {
09161
09162
09163 return (TH1*)gDirectory->Get(hname);
09164 }