loon: UtilHist.cxx Source File

00001 #include "TMath.h"
00002 #include "TH1.h"
00003 #include "TFolder.h"
00004 #include "TClass.h"
00005 
00006 #include "HistMan/HistMan.h"
00007 #include "MessageService/MsgService.h"
00008 #include "Util/UtilMath.h"
00009 
00010 #include "AtNuUtils/UtilHist.h"
00011 
00012 CVSID("$Id: UtilHist.cxx,v 1.3 2007/01/15 19:52:00 rhatcher Exp $");
00013 
00014 //Return a cloned histogram with the figure of merit in each bin as a
00015 //function of a low cut made on the left edge of that bin
00016 TH1* UtilHist::FOMCutLow(TH1 *Signal, TH1 *Background, int iter) {
00017   if(UtilHist::HistMatch(Signal, Background) != 0) return 0;
00018 
00019   Int_t NBins = Signal->GetNbinsX();
00020 
00021   TH1 *FOMOut = dynamic_cast<TH1*>
00022     (Signal->Clone(Form("FOMCutLow%d%s", iter, Signal->GetName())));
00023   FOMOut->SetLineColor(kBlack);
00024   FOMOut->GetYaxis()->SetTitle("FOM Low Cut");
00025   FOMOut->SetBit(TH1::kCanRebin, false);
00026   FOMOut->Reset("ICE");
00027 
00028   //Start with the overflow bin
00029   double NSignal = Signal->GetBinContent(NBins+1);
00030   double NSignalErrSq = Signal->GetBinError(NBins+1) *
00031                         Signal->GetBinError(NBins+1);
00032   double NBackground = Background->GetBinContent(NBins+1);
00033   double NBackgroundErrSq = Background->GetBinError(NBins+1) *
00034                             Background->GetBinError(NBins+1);
00035   for (int i=NBins; i>0; i--) {
00036     NSignal += Signal->GetBinContent(i);
00037     NSignalErrSq += Signal->GetBinError(i) * Signal->GetBinError(i);
00038     NBackground += Background->GetBinContent(i);
00039     NBackgroundErrSq += Background->GetBinError(i) *
00040                         Background->GetBinError(i);
00041     if(NSignal <= 0) continue;
00042     else if(NBackground <= 0) {
00043       FOMOut->SetBinContent(i, 0.);
00044       FOMOut->SetBinError(i, 0.);
00045     }
00046     else {
00047       double StoBG = NSignal / NBackground;
00048       FOMOut->SetBinContent(i, NSignal*StoBG);
00049       FOMOut->SetBinError(i, TMath::Sqrt(StoBG*StoBG*
00050         (4*NSignalErrSq + StoBG*StoBG*NBackgroundErrSq)));
00051     }
00052   }
00053   
00054   return(FOMOut);
00055 }
00056 
00057 //Return a cloned histogram with the figure of merit in each bin as a
00058 //function of a high cut made on the left edge of that bin
00059 TH1* UtilHist::FOMCutHigh(TH1 *Signal, TH1 *Background, int iter) {
00060   if(UtilHist::HistMatch(Signal, Background) != 0) return 0;
00061 
00062   TH1 *FOMOut = dynamic_cast<TH1*>
00063     (Signal->Clone(Form("FOMCutHigh%d%s", iter, Signal->GetName())));
00064   FOMOut->SetLineColor(kBlack);
00065   FOMOut->GetYaxis()->SetTitle("FOM High Cut");
00066   FOMOut->SetBit(TH1::kCanRebin, false);
00067   FOMOut->Reset("ICE");
00068 
00069   //Start with the underflow bin
00070   double NSignal = Signal->GetBinContent(0);
00071   double NSignalErrSq = Signal->GetBinError(0) *
00072                         Signal->GetBinError(0);
00073   double NBackground = Background->GetBinContent(0);
00074   double NBackgroundErrSq = Background->GetBinError(0) *
00075                             Background->GetBinError(0);
00076   for (int i=1; i<Signal->GetNbinsX()+1; i++) {
00077     NSignal += Signal->GetBinContent(i);
00078     NSignalErrSq += Signal->GetBinError(i) * Signal->GetBinError(i);
00079     NBackground += Background->GetBinContent(i);
00080     NBackgroundErrSq += Background->GetBinError(i) *
00081                         Background->GetBinError(i);
00082     if(NSignal <= 0) continue;
00083     else if(NBackground <= 0) {
00084       FOMOut->SetBinContent(i, -0.);
00085       FOMOut->SetBinError(i, -0.);
00086     }
00087     else {
00088       double StoBG = NSignal / NBackground;
00089       FOMOut->SetBinContent(i, NSignal*StoBG);
00090       FOMOut->SetBinError(i, TMath::Sqrt(StoBG*StoBG*
00091         (4*NSignalErrSq + StoBG*StoBG*NBackgroundErrSq)));
00092     }
00093   }
00094   
00095   return(FOMOut);
00096 }
00097 
00098 //Return a cloned histogram with the scaled figure of merit in each bin
00099 //as a function of a low cut made on the left edge of that bin
00100 TH1* UtilHist::SFOMCutLow(TH1 *Signal, TH1 *Background,
00101                           double SignalScale,
00102                           double BackgroundScale, int iter) {
00103   if(UtilHist::HistMatch(Signal, Background) != 0) return 0;
00104 
00105   Int_t NBins = Signal->GetNbinsX();
00106 
00107   TH1 *FOMOut = dynamic_cast<TH1*>
00108     (Signal->Clone(Form("SFOMCutLow%d%s", iter, Signal->GetName())));
00109   FOMOut->SetLineColor(kBlack);
00110   FOMOut->GetYaxis()->SetTitle("Scaled FOM Low Cut");
00111   FOMOut->SetBit(TH1::kCanRebin, false);
00112   FOMOut->Reset("ICE");
00113 
00114   //Start with the overflow bin
00115   double NSignal = SignalScale * Signal->GetBinContent(NBins+1);
00116   double NSignalErrSq = Signal->GetBinError(NBins+1) * SignalScale *
00117                         Signal->GetBinError(NBins+1) * SignalScale;
00118   double NBackground = BackgroundScale * Background->GetBinContent(NBins+1);
00119   double NBackgroundErrSq = Background->GetBinError(NBins+1) * BackgroundScale *
00120                             Background->GetBinError(NBins+1) * BackgroundScale;
00121   for (int i=NBins; i>0; i--) {
00122     NSignal += SignalScale * Signal->GetBinContent(i);
00123     NSignalErrSq += SignalScale * Signal->GetBinError(i) *
00124                     SignalScale * Signal->GetBinError(i);
00125     NBackground += BackgroundScale * Background->GetBinContent(i);
00126     NBackgroundErrSq += Background->GetBinError(i) * BackgroundScale *
00127                         Background->GetBinError(i) * BackgroundScale; 
00128     if(NSignal <= 0) continue;
00129     else if (NBackground <= 0) {
00130       FOMOut->SetBinContent(i, -0.);
00131       FOMOut->SetBinError(i, -0.);
00132     }
00133     else {
00134       FOMOut->SetBinContent(i, NSignal*NSignal/(NSignal+NBackground));
00135       FOMOut->SetBinError(i, TMath::Sqrt(
00136          ((4*NSignal*NSignal*NBackground /
00137            ((NSignal+NBackground) * (NSignal+NBackground) *
00138             (NSignal+NBackground))) * NSignalErrSq) +
00139          ((4*NSignal*NSignal*NSignal*NSignal /
00140            ((NSignal+NBackground) * (NSignal+NBackground) *
00141             (NSignal+NBackground) * (NSignal+NBackground))
00142             ) * (NSignalErrSq+NBackgroundErrSq))
00143         ));
00144     }
00145   }
00146   
00147   
00148   return(FOMOut);
00149 }
00150 
00151 //Return a cloned histogram with the scaled figure of merit in each bin
00152 //as a function of a high cut made on the left edge of that bin
00153 TH1* UtilHist::SFOMCutHigh(TH1 *Signal, TH1 *Background,
00154                            double SignalScale,
00155                            double BackgroundScale, int iter) {
00156   if(UtilHist::HistMatch(Signal, Background) != 0) return 0;
00157 
00158   TH1 *FOMOut = dynamic_cast<TH1*>
00159     (Signal->Clone(Form("SFOMCutHigh%d%s", iter, Signal->GetName())));
00160   FOMOut->SetLineColor(kBlack);
00161   FOMOut->GetYaxis()->SetTitle("Scaled FOM High Cut");
00162   FOMOut->SetBit(TH1::kCanRebin, false);
00163   FOMOut->Reset("ICE");
00164 
00165   //Start with the underflow bin
00166   double NSignal = SignalScale * Signal->GetBinContent(0);
00167   double NSignalErrSq = Signal->GetBinError(0) * SignalScale *
00168                         Signal->GetBinError(0) * SignalScale;
00169   double NBackground = BackgroundScale * Background->GetBinContent(0);
00170   double NBackgroundErrSq = Background->GetBinError(0) * BackgroundScale *
00171                             Background->GetBinError(0) * BackgroundScale;
00172   for (int i=1; i<Signal->GetNbinsX()+1; i++) {
00173     NSignal += SignalScale * Signal->GetBinContent(i);
00174     NSignalErrSq += SignalScale * Signal->GetBinError(i) *
00175                     SignalScale * Signal->GetBinError(i);
00176     NBackground += BackgroundScale * Background->GetBinContent(i);
00177     NBackgroundErrSq += Background->GetBinError(i) * BackgroundScale *
00178                         Background->GetBinError(i) * BackgroundScale; 
00179     if(NSignal <= 0) continue;
00180     else if (NBackground <= 0) {
00181       FOMOut->SetBinContent(i, -0.);
00182       FOMOut->SetBinError(i, -0.);
00183     }
00184     else {
00185       FOMOut->SetBinContent(i, NSignal*NSignal/(NSignal+NBackground));
00186       FOMOut->SetBinError(i, TMath::Sqrt(
00187          ((4*NSignal*NSignal*NBackground /
00188            ((NSignal+NBackground) * (NSignal+NBackground) *
00189             (NSignal+NBackground))) * NSignalErrSq) +
00190          ((4*NSignal*NSignal*NSignal*NSignal /
00191            ((NSignal+NBackground) * (NSignal+NBackground) *
00192             (NSignal+NBackground) * (NSignal+NBackground))
00193             ) * (NSignalErrSq+NBackgroundErrSq))
00194         ));
00195 /*
00196           (NSignal*NSignal/(NSignal+NBackground)) * (
00197             (4*NSignalErrSq/(NSignal*NSignal)) +
00198             ((NSignalErrSq+NBackgroundErrSq)/
00199               ((NSignal+NBackground)*(NSignal+NBackground)))
00200             )
00201           ));
00202 */
00203     }
00204   }
00205   
00206   return(FOMOut);
00207 }
00208 
00209 void UtilHist::ScaleHists(Float_t AtmScale, Float_t CosScale,
00210                           TFolder *tfol) {
00211   //For empty parameter, start at the top
00212   TFolder *tf = tfol;
00213   if (!tf) {
00214     HistMan hm("");
00215     tf = &(hm.BaseFolder());
00216   }
00217 
00218   string CosSuf = "Cosmic";
00219   string AtmSuf = "Atmos";
00220 
00221   TCollection *flist = tf->GetListOfFolders();
00222   TIter itr(flist->MakeIterator());
00223   TObject *fobj = 0;
00224   while( (fobj=itr()) ) {
00225     if (fobj->InheritsFrom("TFolder")) {
00226       ScaleHists(AtmScale,CosScale,dynamic_cast<TFolder*>(fobj));
00227       continue;
00228     }
00229     if (fobj->InheritsFrom("TH1")) {
00230       TH1 *thishist = dynamic_cast<TH1*>(fobj);
00231       string hname = thishist->GetName();
00232       if (hname.length() > CosSuf.length() &&
00233           hname.find(CosSuf) == hname.length()-CosSuf.length()) {
00234         thishist->Scale(CosScale);
00235       }
00236 
00237       if (hname.length() > AtmSuf.length() &&
00238           hname.find(AtmSuf) == hname.length()-AtmSuf.length()) {
00239         thishist->Scale(AtmScale);
00240       }
00241     }
00242   }
00243 
00244   return;
00245 }
00246 
00247 int UtilHist::HistMatch(TH1 *H1, TH1 *H2) {
00248   // Make sure two histograms match for comparison
00249   // writes detailed warnings to cout if there's a problem.
00250   // Status flag:
00251   //   0 = OK
00252   //   1 = One or both of the histograms do not exist
00253   //   2 = Different histogram types
00254   //   3 = Different number of bins
00255   //   4 = Different axes ranges
00256   //   5 = One of both of the histograms has no entries
00257   //
00258 
00259   //One or both histograms do not exist
00260   if(!H1 || !H2) return(1);
00261 
00262   TClass* cl1 = H1->Class();
00263   TClass* cl2 = H2->Class();
00264   //Different types of histogram
00265   if(!(cl1->InheritsFrom(cl2)) || !(cl2->InheritsFrom(cl1))) return(2);
00266 
00267   //Different number of bins in each histogram
00268   if((H1->GetNbinsX() != H2->GetNbinsX()) ||
00269      (H1->GetNbinsY()!=H2->GetNbinsY()) ||
00270      (H1->GetNbinsZ()!=H2->GetNbinsZ()) ) return(3);
00271 
00272   //Histogram range mis-match
00273   if((H1->GetXaxis()->GetXmin() != H2->GetXaxis()->GetXmin()) ||
00274      (H1->GetXaxis()->GetXmax() != H2->GetXaxis()->GetXmax()) ||
00275      (H1->GetYaxis()->GetXmin() != H2->GetYaxis()->GetXmin()) ||
00276      (H1->GetYaxis()->GetXmax() != H2->GetYaxis()->GetXmax()) ||
00277      (H1->GetZaxis()->GetXmin() != H2->GetZaxis()->GetXmin()) ||
00278      (H1->GetZaxis()->GetXmax() != H2->GetZaxis()->GetXmax()))
00279     return(4);
00280 
00281 
00282   //One or both of the histograms are empty
00283   if(H1->GetEntries() <= 0. || H2->GetEntries() <= 0.) return(5);
00284 
00285   //Checks out
00286   return 0;
00287 }
00288 
00289 //Have a Min and Max function to apply to TH1 objects, since the
00290 //GetMinimum and GetMaximum methods don't always work.
00291 double UtilHist::HMin(const TH1 *hist, double min) {
00292   if(!hist) return(min);
00293   double lmin = DBL_MAX;
00294   for (int i=1; i<=hist->GetNbinsX(); i++) {
00295     double bcon = hist->GetBinContent(i);
00296     if(! UtilMath::sameValue(hist->GetBinContent(i),
00297                              hist->GetBinError(i)) )
00298       bcon -= hist->GetBinError(i);
00299     if(bcon < lmin && bcon > min) lmin = bcon;
00300   }
00301   return lmin;
00302 }
00303 
00304 double UtilHist::HMax(const TH1 *hist, double max) {
00305   if(!hist) return(max);
00306   double lmax = -DBL_MAX;
00307   for (int i=1; i<=hist->GetNbinsX(); i++) {
00308     double bcon = hist->GetBinContent(i);
00309     bcon += hist->GetBinError(i);
00310     if(bcon > lmax && bcon < max) lmax = bcon;
00311   }
00312   return lmax;
00313 }
00314 
00315 double UtilHist::HMin(const vector<TH1*> &hists, double min) {
00316   MSG("Util",Msg::kVerbose) << "HMin on " << hists.size() << " histograms" << endl;
00317   MSG("Util",Msg::kVerbose) << " min = " << min << endl;
00318   if(hists.size() == 0) return(min);
00319   double lmin = DBL_MAX;
00320   for (unsigned int i=0; i<hists.size(); i++) {
00321     MSG("Util",Msg::kVerbose) << " Hist Name = " << hists[0]->GetName() << endl;
00322     double llmin = UtilHist::HMin(hists[i], min);
00323     MSG("Util",Msg::kVerbose) << "  llmin = " << llmin << endl;
00324     if(llmin < lmin) lmin = llmin;
00325   }
00326   MSG("Util",Msg::kVerbose) << "lmin = " << lmin << endl;
00327   return lmin;
00328 }
00329 
00330 double UtilHist::HMax(const vector<TH1*> &hists, double max) {
00331   MSG("Util",Msg::kVerbose) << "HMax on " << hists.size() << " histograms" << endl;
00332   MSG("Util",Msg::kVerbose) << " max = " << max << endl;
00333   if(hists.size() == 0) return(max);
00334   double lmax = -DBL_MAX;
00335   for (unsigned int i=0; i<hists.size(); i++) {
00336     MSG("Util",Msg::kVerbose) << "Hist Name = " << hists[0]->GetName() << endl;
00337     double llmax = UtilHist::HMax(hists[i], max);
00338     MSG("Util",Msg::kVerbose) << "  llmax = " << llmax << endl;
00339     if(llmax > lmax) lmax = llmax;
00340   }
00341   MSG("Util",Msg::kVerbose) << "lmax = " << lmax << endl;
00342   return lmax;
00343 }