SFH: RegH1F.C Source File

00001 
00030 #include "jbltools/sfh/RegH1F.h"
00031 #include "jbltools/sfh/Binning.h"
00032 
00033 // system headers
00034 #include <iostream>
00035 #include <cassert>
00036 
00037 // root headers
00038 #include <TFile.h>
00039 #include <TF1.h>
00040 
00041 
00042 using std::endl;
00043 using std::cout;
00044 using std::cerr;
00045 
00046 static const char *ident="@(#)$Id: RegH1F.C,v 1.15 2005/11/11 17:21:16 rweber Exp $";
00047 
00048 
00049 RegH1F::RegH1F (const char* name, 
00050                 const char* title, 
00051                 Int_t nbinsx, Axis_t xlow, Axis_t xup,
00052                 const ROListPoR& hhl) 
00053 : TH1F (name, title, nbinsx, xlow, xup), 
00054   RegO (hhl)
00055 {
00056   if(fSumw2.fN == 0) Sumw2();
00057 }
00058   
00059   
00060 RegH1F::RegH1F (const char* name, 
00061                 const char* title, 
00062                 Int_t nbinsx, const Float_t* xbins,
00063                 const ROListPoR& hhl) 
00064 : TH1F (name, title, nbinsx, xbins), 
00065   RegO (hhl)
00066 {
00067   if(fSumw2.fN == 0) Sumw2();
00068 }
00069 
00070 RegH1F::RegH1F (const char* name, 
00071                 const char* title, 
00072                 Int_t nbinsx, const Double_t* xbins,
00073                 const ROListPoR& hhl) 
00074 : TH1F (name, title, nbinsx, xbins), 
00075   RegO (hhl)  
00076 {
00077   if(fSumw2.fN == 0) Sumw2();
00078 }
00079 
00080 RegH1F::RegH1F (const char* name, 
00081                 const char* title, 
00082                 const Binning& binning,
00083                 const ROListPoR& hhl) 
00084 : TH1F (name, title, binning.getNBins(), binning.getLowerEdge(), binning.getUpperEdge()), 
00085   RegO (hhl)  
00086 {
00087   // For an equidistant binning, everything is OK already;
00088   // for a non-equidistant binning, the number of bins is fine,
00089   // but we have to adjust the bin edges
00090   if (!binning.isEquidistant()) {
00091     TAxis *xaxis = GetXaxis();
00092     assert (xaxis);
00093     assert (binning.getEdges());
00094     xaxis->Set (binning.getNBins(), binning.getEdges());
00095   }
00096   if (binning.hasBinLabels()) {
00097     TAxis *xaxis = GetXaxis();
00098     assert (xaxis);
00099     for (int i = 0; i< binning.getNBins(); ++i) {
00100       if (const char *label = binning.getBinLabel (i)) {
00101         xaxis->SetBinLabel (i+1, label);
00102       }
00103     }
00104     LabelsOption ("u", "X");
00105   }
00106   if(fSumw2.fN == 0) Sumw2();
00107 }
00108 
00109 
00110 RegH1F::RegH1F (const char* name, TFile& file, const ROListPoR& hhl) 
00111 : RegO (hhl)
00112 {
00113   TH1F *h = dynamic_cast<TH1F *>(file.Get(name));
00114   if (h) {
00115     h->Copy(*this);
00116     cout << "RegH1F::RegH1F (const char* name, TFile& file): "
00117          << "histogram " << name << " read from " << file.GetName() << "!" << endl;
00118   }
00119   else {
00120     cerr << "RegH1F::RegH1F (const char* name, TFile& file): "
00121          << "histogram " << name << " not found in file " << file.GetName() << "!" << endl;
00122   }
00123   delete h;
00124   if(fSumw2.fN == 0) Sumw2();
00125 }
00126 
00127 RegH1F::RegH1F (const char* name, const char* title, 
00128                 const TH1F& source, const ROListPoR& hhl) 
00129 : TH1F (source), RegO (hhl) {
00130   if (name) SetName (name);
00131   if (title) SetTitle (title);
00132   if(fSumw2.fN == 0) Sumw2();
00133 }
00134 
00135 RegH1F::~RegH1F() {}
00136 
00137 
00138 Stat_t RegH1F::GetSumOfErrors2() const
00139 {
00140 //   -*-*-*-*-*-*Return the sum of error squares excluding under/overflows*-*-*-*-*
00141 //               ===================================================
00142   Int_t bin,binx,biny,binz;
00143   Stat_t sum =0;
00144   for(binz=1; binz<=fZaxis.GetNbins(); binz++) {
00145      for(biny=1; biny<=fYaxis.GetNbins(); biny++) {
00146         for(binx=1; binx<=fXaxis.GetNbins(); binx++) {
00147            bin = GetBin(binx,biny,binz);
00148 //           sum += GetBinError2(bin);
00149            sum += GetBinError(bin)*GetBinError(bin);
00150         }
00151      }
00152   }
00153   return sum;
00154 }
00155 void RegH1F::getHistInfo (Option_t* option,   
00156                   Stat_t& content,    
00157                   Stat_t& error2      
00158                   ) const {
00159   content = 0;
00160   error2 = 0;
00161   
00162   TString opt = option;
00163   opt.ToLower();
00164   
00165   if (opt.Contains("m")) {
00166     content = GetMean();
00167     error2 = GetRMS();
00168     error2 *= error2;
00169   }
00170   else {
00171     content = GetSumOfWeights();
00172     error2 = GetSumOfErrors2();
00173   }
00174 }            
00175 
00176 
00177 void RegH1F::Divide (TF1* f1, Double_t c1) {
00178   TH1F::Divide (f1, c1);
00179 }
00180 
00181 void RegH1F::Divide (const TH1* h1) {
00182   TH1F::Divide (h1);
00183 }
00184 
00185 void RegH1F::Divide(const TH1 *h1, const TH1 *h2, Double_t c1, Double_t c2, Option_t *option) {
00186 // THIS is the original TH1::Divide method, 
00187 // but errors are done correctly for weighted events!
00188 // (and c1 and c2 are ignored for 'B' option)
00189 
00190 // BL 13.6.03: Option "s" signals that h1 contains a (weighted) subsample
00191 // of events of h2, errors are then calculated accordingly
00192 // For this option, it is assumed that all events in h2 originally had the same
00193 // weight W=e2/b2, so the number of events in one bin of h2 is assumed to be b2*b2/e2.
00194 
00195    TString opt = option;
00196    opt.ToLower();
00197    Bool_t binomial = kFALSE;
00198    if (opt.Contains("b")) binomial = kTRUE;
00199    Bool_t subsample = kFALSE;
00200    if (opt.Contains("s")) {
00201      subsample = kTRUE;
00202      // std::cerr << "RegH1F::Divide : using option 's'\n";
00203    }
00204    if (!h1 || !h2) {
00205       Error("Divide","Attempt to divide by a non-existing histogram");
00206       return;
00207    }
00208 
00209    Int_t nbinsx = GetNbinsX();
00210    Int_t nbinsy = GetNbinsY();
00211    Int_t nbinsz = GetNbinsZ();
00212 //   - Check histogram compatibility
00213    if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
00214     || nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()) {
00215       Error("Divide","Attempt to divide histograms with different number of bins");
00216       return;
00217    }
00218    if (!c2) {
00219       Error("Divide","Coefficient of dividing histogram cannot be zero");
00220       return;
00221    }
00222 //   - Issue a Warning if histogram limits are different
00223    if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
00224        GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
00225        GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
00226        GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
00227        GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
00228        GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
00229        Warning("Divide","Attempt to divide histograms with different axis limits");
00230    }
00231    if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
00232        GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
00233        GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
00234        GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
00235        GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
00236        GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
00237        Warning("Divide","Attempt to divide histograms with different axis limits");
00238    }
00239    if (fDimension < 2) nbinsy = -1;
00240    if (fDimension < 3) nbinsz = -1;
00241 
00242 //    Create Sumw2 if h1 or h2 have Sumw2 set
00243    if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
00244 
00245 //   - Reset statistics
00246    fEntries = fTsumw = 0;
00247 
00248 //   - Loop on bins (including underflows/overflows)
00249    Int_t bin, binx, biny, binz;
00250    Double_t b1,b2,w,d1,d2;
00251    d1 = c1*c1;
00252    d2 = c2*c2;
00253 //    std::cout << "RegisteredHisto::Divide: treating histogram " << GetTitle() << endl;
00254    for (binz=0;binz<=nbinsz+1;binz++) {
00255       for (biny=0;biny<=nbinsy+1;biny++) {
00256          for (binx=0;binx<=nbinsx+1;binx++) {
00257             bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
00258             b1  = h1->GetBinContent(bin);
00259             b2  = h2->GetBinContent(bin);
00260             if (b2) {
00261               if (binomial) w = b1/b2;
00262               else w = c1*b1/(c2*b2);
00263             } 
00264             else  {w = 0;}
00265             SetBinContent(bin,w);
00266             fEntries++;
00267             if (fSumw2.fN) {
00268                if (!b2) { fSumw2.fArray[bin] = 0; continue;}
00269                Double_t e1 = h1->GetBinError(bin);
00270                Double_t e2 = h2->GetBinError(bin);
00271                if (binomial) {
00272                   //fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/(c2*b2));//this is the formula in Hbook/Hoper1
00273                   fSumw2.fArray[bin] = TMath::Abs(w*(1-w)/b2);
00274                } else if (subsample) {
00275                   fSumw2.fArray[bin] =  d1*(e1*e1*b2*(b2 - 2*b1) + e2*e2*b1*b1)/(d2*b2*b2*b2*b2);
00276                } else {
00277                   Double_t b22= b2*b2*d2;
00278                   fSumw2.fArray[bin] = d1*d2*(e1*e1*b2*b2 + e2*e2*b1*b1)/(b22*b22);
00279                }
00280             }
00281          }
00282       }
00283    }
00284    Stat_t s[10];
00285    GetStats(s);
00286    PutStats(s);
00287 }
00288 
00289 void RegH1F::Divide(const TH1 *h1, const TH1 *h2, const TH1 *h3, Double_t c1, Double_t c2) {
00290 
00291    if (!h1 || !h2 || !h3) {
00292       Error("Divide","Attempt to divide by a non-existing histogram");
00293       return;
00294    }
00295 
00296    Int_t nbinsx = GetNbinsX();
00297    Int_t nbinsy = GetNbinsY();
00298    Int_t nbinsz = GetNbinsZ();
00299 //   - Check histogram compatibility
00300    if (nbinsx != h1->GetNbinsX() || nbinsy != h1->GetNbinsY() || nbinsz != h1->GetNbinsZ()
00301     || nbinsx != h2->GetNbinsX() || nbinsy != h2->GetNbinsY() || nbinsz != h2->GetNbinsZ()
00302     || nbinsx != h3->GetNbinsX() || nbinsy != h3->GetNbinsY() || nbinsz != h3->GetNbinsZ()) {
00303       Error("Divide","Attempt to divide histograms with different number of bins");
00304       return;
00305    }
00306    if (!c2) {
00307       Error("Divide","Coefficient of dividing histogram cannot be zero");
00308       return;
00309    }
00310 //   - Issue a Warning if histogram limits are different
00311    if (GetXaxis()->GetXmin() != h1->GetXaxis()->GetXmin() ||
00312        GetXaxis()->GetXmax() != h1->GetXaxis()->GetXmax() ||
00313        GetYaxis()->GetXmin() != h1->GetYaxis()->GetXmin() ||
00314        GetYaxis()->GetXmax() != h1->GetYaxis()->GetXmax() ||
00315        GetZaxis()->GetXmin() != h1->GetZaxis()->GetXmin() ||
00316        GetZaxis()->GetXmax() != h1->GetZaxis()->GetXmax()) {
00317        Warning("Divide","Attempt to divide histograms with different axis limits");
00318    }
00319    if (GetXaxis()->GetXmin() != h2->GetXaxis()->GetXmin() ||
00320        GetXaxis()->GetXmax() != h2->GetXaxis()->GetXmax() ||
00321        GetYaxis()->GetXmin() != h2->GetYaxis()->GetXmin() ||
00322        GetYaxis()->GetXmax() != h2->GetYaxis()->GetXmax() ||
00323        GetZaxis()->GetXmin() != h2->GetZaxis()->GetXmin() ||
00324        GetZaxis()->GetXmax() != h2->GetZaxis()->GetXmax()) {
00325        Warning("Divide","Attempt to divide histograms with different axis limits");
00326    }
00327    if (GetXaxis()->GetXmin() != h3->GetXaxis()->GetXmin() ||
00328        GetXaxis()->GetXmax() != h3->GetXaxis()->GetXmax() ||
00329        GetYaxis()->GetXmin() != h3->GetYaxis()->GetXmin() ||
00330        GetYaxis()->GetXmax() != h3->GetYaxis()->GetXmax() ||
00331        GetZaxis()->GetXmin() != h3->GetZaxis()->GetXmin() ||
00332        GetZaxis()->GetXmax() != h3->GetZaxis()->GetXmax()) {
00333        Warning("Divide","Attempt to divide histograms with different axis limits");
00334    }
00335    if (fDimension < 2) nbinsy = -1;
00336    if (fDimension < 3) nbinsz = -1;
00337 
00338 //    Create Sumw2 if h1 or h2 have Sumw2 set
00339    if (fSumw2.fN == 0 && (h1->GetSumw2N() != 0 || h2->GetSumw2N() != 0)) Sumw2();
00340 
00341 //   - Reset statistics
00342    fEntries = fTsumw = 0;
00343 
00344 //   - Loop on bins (including underflows/overflows)
00345    Int_t bin, binx, biny, binz;
00346    Double_t b1,b2,w,d1,d2;
00347    d1 = c1*c1;
00348    d2 = c2*c2;
00349 //    std::cout << "RegisteredHisto::Divide: treating histogram " << GetTitle() << endl;
00350    for (binz=0;binz<=nbinsz+1;binz++) {
00351       for (biny=0;biny<=nbinsy+1;biny++) {
00352          for (binx=0;binx<=nbinsx+1;binx++) {
00353             bin = binx +(nbinsx+2)*(biny + (nbinsy+2)*binz);
00354             b1  = h1->GetBinContent(bin);
00355             b2  = h2->GetBinContent(bin);
00356             if (b2) {
00357               w = c1*b1/(c2*b2);
00358             } 
00359             else  {w = 0;}
00360             SetBinContent(bin,w);
00361             fEntries++;
00362             if (fSumw2.fN) {
00363                if (!b2) { fSumw2.fArray[bin] = 0; continue;}
00364                Double_t e1 = h1->GetBinError(bin);
00365                Double_t e2 = h2->GetBinError(bin);
00366                Double_t b3 = h3->GetBinContent(bin);
00367                
00368                fSumw2.fArray[bin] = d1*(e1*e1*b2*b2 - 2*b3*b1*b2 + e2*e2*b1*b1)/(d2*b2*b2*b2*b2);
00369             }
00370          }
00371       }
00372    }
00373    Stat_t s[10];
00374    GetStats(s);
00375    PutStats(s);
00376 }
00377 
00378 // DrawFunc by Ronnie Weber
00379 void RegH1F::DrawFunc (Option_t* option, const char* fname) {
00380   if(TF1 *f = this->GetFunction(fname)) {
00381     std::cout << "RegH1F::DrawFunc: temp. enabling function " << fname 
00382               << " for drawing of " << this->GetName() << "\n";
00383     f->SetBit(BIT(9), 0);
00384   } else {
00385     std::cout << "RegH1F::DrawFunc: ERROR - could not find function " << fname
00386               << " for drawing of " << this->GetName() << "\n";
00387   }
00388   this->DrawCopy(option);
00389   if(TF1 *f = this->GetFunction(fname)) {
00390 //     std::cout << "DrawFunc: Disabling function " << fname 
00391 //               << " for drawing of " << this->GetName() << ".\n";
00392     f->SetBit(BIT(9), 1);
00393   }
00394 }
00395