| 1 |
#ifndef GENERIC_PLOTTER_H
|
| 2 |
#define GENERIC_PLOTTER_H
|
| 3 |
|
| 4 |
#include <stdio.h>
|
| 5 |
#include <iostream>
|
| 6 |
#include <vector>
|
| 7 |
#include <string>
|
| 8 |
#include <TLegend.h>
|
| 9 |
#include <TAxis.h>
|
| 10 |
#include <TH1.h>
|
| 11 |
#include <math.h>
|
| 12 |
|
| 13 |
struct coolRatio {
|
| 14 |
//coolRatio() :
|
| 15 |
|
| 16 |
// User must clean up the created ratio histograms.
|
| 17 |
//TH1 make_rebinned_ratios(TH1* theHist, TH1* theReference, double maxUncertainty);
|
| 18 |
//static double ratioError2(double numerator, double numeratorError2, double denominator, double denominatorError2);
|
| 19 |
|
| 20 |
|
| 21 |
|
| 22 |
TH1* make_rebinned_ratios(TH1* theHist, TH1* theReference, double maxUncertainty, bool useReferenceError = false, const int number = 0)
|
| 23 |
{
|
| 24 |
std::vector<TH1*> hist;
|
| 25 |
std::vector<TH1*> ratios;
|
| 26 |
hist.push_back(theHist);
|
| 27 |
hist.push_back(theReference);
|
| 28 |
const int reference = 1;
|
| 29 |
|
| 30 |
const int numHistograms = int(hist.size());
|
| 31 |
assert(numHistograms >= 2);
|
| 32 |
ratios.clear();
|
| 33 |
|
| 34 |
// Find nonzero range
|
| 35 |
const TH1* denominator = hist[reference];
|
| 36 |
const TAxis* rebinAxis = denominator->GetXaxis();
|
| 37 |
const int binBound = rebinAxis->GetNbins() + 1; // overflow
|
| 38 |
int firstBin = 0;
|
| 39 |
for (bool allEmpty = true; allEmpty && ++firstBin < binBound;)
|
| 40 |
for (int i = 0; allEmpty && i < numHistograms; ++i)
|
| 41 |
allEmpty &= (hist[i] == 0 || hist[i]->GetBinContent(firstBin) == 0);
|
| 42 |
int lastBin = binBound;
|
| 43 |
for (bool allEmpty = true; allEmpty && --lastBin > 0;)
|
| 44 |
for (int i = 0; allEmpty && i < numHistograms; ++i)
|
| 45 |
allEmpty &= (hist[i] == 0 || hist[i]->GetBinContent(lastBin) == 0);
|
| 46 |
if (lastBin <= firstBin) return theHist; // no content
|
| 47 |
|
| 48 |
|
| 49 |
// Partition bins in range
|
| 50 |
std::vector<double> binEdges;
|
| 51 |
if (firstBin > 1) binEdges.push_back(rebinAxis->GetBinLowEdge(1));
|
| 52 |
double previousEdge = rebinAxis->GetBinLowEdge(firstBin);
|
| 53 |
std::vector<double> zeros (numHistograms,0);
|
| 54 |
std::vector<std::vector<double> > sumNum (1, zeros);
|
| 55 |
std::vector<std::vector<double> > sumNumErr2(1, zeros);
|
| 56 |
std::vector<double> sumDen (1, 0);
|
| 57 |
std::vector<double> sumDenErr2(1, 0);
|
| 58 |
|
| 59 |
////std::cout << std::endl << rebinAxis->GetXmin() << " ( " << rebinAxis->GetBinLowEdge(firstBin) << " , ";
|
| 60 |
////std::cout << rebinAxis->GetBinUpEdge(lastBin) << " ) " << rebinAxis->GetXmax() << std::endl;
|
| 61 |
for (int bin = firstBin; bin <= lastBin; ++bin) {
|
| 62 |
// Find a filled block
|
| 63 |
sumDen .back() += denominator->GetBinContent(bin);
|
| 64 |
sumDenErr2.back() += pow(denominator->GetBinError(bin), 2);
|
| 65 |
bool aboveThreshold = (sumDen.back() != 0); // denominator must not vanish
|
| 66 |
|
| 67 |
for (int iHist = 0; iHist < numHistograms; ++iHist) {
|
| 68 |
if (iHist == reference) continue;
|
| 69 |
if (hist[iHist] == 0) continue;
|
| 70 |
double num = hist[iHist]->GetBinContent(bin);
|
| 71 |
double numErr2 = pow(hist[iHist]->GetBinError(bin), 2);
|
| 72 |
sumNum .back()[iHist] += num;
|
| 73 |
sumNumErr2.back()[iHist] += numErr2;
|
| 74 |
|
| 75 |
const double error2 = ratioError2(sumNum.back()[iHist],sumNumErr2.back()[iHist],sumDen.back(),sumDenErr2.back(),useReferenceError);
|
| 76 |
aboveThreshold &= (sqrt(error2) * sumDen.back() < maxUncertainty * sumNum.back()[iHist]);
|
| 77 |
if (maxUncertainty > 100.) aboveThreshold = true;
|
| 78 |
} // end loop over histograms
|
| 79 |
|
| 80 |
if (aboveThreshold && bin < lastBin) {
|
| 81 |
//std::cout << " + [" << previousEdge << " , " << rebinAxis->GetBinUpEdge(bin) << "] = " << sumDen.back() << std::endl;
|
| 82 |
binEdges.push_back(previousEdge);
|
| 83 |
sumNum .push_back(zeros); sumNumErr2.push_back(zeros);
|
| 84 |
sumDen .push_back(0); sumDenErr2.push_back(0);
|
| 85 |
previousEdge = rebinAxis->GetBinUpEdge(bin);
|
| 86 |
}
|
| 87 |
} // end loop over bins in reference
|
| 88 |
|
| 89 |
if (binEdges.empty() || binEdges.back()!=previousEdge) binEdges.push_back(previousEdge);
|
| 90 |
binEdges.push_back(rebinAxis->GetBinUpEdge (lastBin ));
|
| 91 |
if (lastBin < binBound - 1) binEdges.push_back(rebinAxis->GetBinLowEdge(binBound));
|
| 92 |
|
| 93 |
////for (unsigned i=0; i<binEdges.size(); ++i) std::cout << " " << binEdges[i]; std::cout << std::endl;
|
| 94 |
////for (unsigned i=0; i<sumDen.size(); ++i) std::cout << " " << sumDen[i]; std::cout << std::endl;
|
| 95 |
|
| 96 |
|
| 97 |
// Make ratios according to new bin edges
|
| 98 |
ratios.reserve(numHistograms - 1);
|
| 99 |
const unsigned numBins = binEdges.size() - 1u;
|
| 100 |
const unsigned numContent = sumNum.size();
|
| 101 |
assert(numBins >= numContent);
|
| 102 |
|
| 103 |
|
| 104 |
for (int iHist = 0; iHist < numHistograms; ++iHist) {
|
| 105 |
if (iHist == reference) continue;
|
| 106 |
if (hist[iHist] == 0) continue;
|
| 107 |
TH1* ratio = static_cast<TH1*>(hist[iHist]->Clone("Ratio"));
|
| 108 |
TH1* refError = static_cast<TH1*>(hist[iHist]->Clone("RefError"));
|
| 109 |
ratio->Reset ();
|
| 110 |
ratio->SetBins (numBins, &(binEdges[0]));
|
| 111 |
refError->Reset ();
|
| 112 |
refError->SetBins (numBins, &(binEdges[0]));
|
| 113 |
if (ratio->GetSumw2N() == 0) ratio->Sumw2();
|
| 114 |
if (refError->GetSumw2N() == 0) refError->Sumw2();
|
| 115 |
TArrayD* error2 = ratio->GetSumw2();
|
| 116 |
TArrayD* refError2 = refError->GetSumw2();
|
| 117 |
for (unsigned index = 0, bin = 1 + (firstBin > 1); index < numContent; ++index, ++bin) {
|
| 118 |
if (sumDen[index] == 0) continue;
|
| 119 |
ratio ->SetBinContent(bin, sumNum[index][iHist] / sumDen[index]);
|
| 120 |
error2->SetAt(ratioError2(sumNum[index][iHist],sumNumErr2[index][iHist],sumDen[index],sumDenErr2[index],useReferenceError), bin);
|
| 121 |
if (sumNum[index][iHist] == 0 && sumDen[index] == 0){
|
| 122 |
refError ->SetBinContent(bin, 0.);
|
| 123 |
refError2->SetAt(0., bin);
|
| 124 |
}
|
| 125 |
else{
|
| 126 |
refError ->SetBinContent(bin, 1.);
|
| 127 |
refError2->SetAt(ratio1Error2(sumDen[index],sumDenErr2[index]), bin);
|
| 128 |
}
|
| 129 |
} // end loop over bins
|
| 130 |
ratios.push_back(ratio);
|
| 131 |
ratios.push_back(refError);
|
| 132 |
} // end loop over histograms
|
| 133 |
return ratios[number];
|
| 134 |
}
|
| 135 |
|
| 136 |
private:
|
| 137 |
double ratio1Error2(double denominator, double denominatorError2)
|
| 138 |
{
|
| 139 |
if (denominator == 0) return 0;
|
| 140 |
return (denominatorError2/denominator/denominator );
|
| 141 |
}
|
| 142 |
double ratioError2(double numerator, double numeratorError2, double denominator, double denominatorError2,bool useDenominatorError=true)
|
| 143 |
{
|
| 144 |
if (numerator == 0) return 0;
|
| 145 |
if (denominator == 0) return 0;
|
| 146 |
if (!useDenominatorError) denominatorError2 = 0;
|
| 147 |
const double ratio = numerator / denominator;
|
| 148 |
return ratio*ratio*( numeratorError2/numerator/numerator + denominatorError2/denominator/denominator );
|
| 149 |
}
|
| 150 |
|
| 151 |
};
|
| 152 |
|
| 153 |
#endif
|
