| 3 |
|
#include <TFile.h> |
| 4 |
|
#include <TH1.h> |
| 5 |
|
#include <TH2.h> |
| 6 |
+ |
#include <THStack.h> |
| 7 |
|
#include <TLatex.h> |
| 8 |
|
#include <TCanvas.h> |
| 9 |
|
#include <TMath.h> |
| 10 |
|
#include <TLegend.h> |
| 11 |
|
#include <TObject.h> |
| 12 |
+ |
#include <TStyle.h> |
| 13 |
|
#include <iostream> |
| 14 |
|
#include <sstream> |
| 15 |
|
#include <fstream> |
| 32 |
|
}; |
| 33 |
|
|
| 34 |
|
//function for doing KS test |
| 35 |
< |
vector<testMC> doKStest(Double_t NsS, Double_t Ns1, Double_t Ns2, TH1F* mixS, TH1F* s1, TH1F* s2) { |
| 35 |
> |
vector<testMC> doKStest(const vector<Double_t>& vnevt, map<string,TH1F*>& mhisto) { |
| 36 |
|
vector<testMC> output; |
| 37 |
|
//define the scale factors |
| 38 |
|
Double_t sf1 = 0.0; // QCD |
| 39 |
|
Double_t sf2 = 0.0; // Wjets |
| 40 |
+ |
Double_t Ns = 0.0; // non-Wjets (ttbar,Zjets,single top...) |
| 41 |
+ |
for (UInt_t i = 3; i < vnevt.size(); ++i) |
| 42 |
+ |
Ns += vnevt[i]; |
| 43 |
+ |
cout << "Total number of Non Wjets events = " << Ns << endl; |
| 44 |
|
//KS test |
| 45 |
|
do { |
| 46 |
< |
sf1 = (NsS - Ns2*sf2)/Ns1; |
| 46 |
> |
sf1 = (vnevt.at(0) - Ns - vnevt.at(2)*sf2)/vnevt.at(1); |
| 47 |
|
if (sf1 < 0) break; |
| 48 |
|
//cout << "..........sf1 = " << sf1 << endl; |
| 49 |
< |
int nbins = mixS->GetNbinsX(); |
| 50 |
< |
double xmin = mixS->GetXaxis()->GetBinLowEdge(1); |
| 51 |
< |
double xmax = mixS->GetXaxis()->GetBinUpEdge(nbins); |
| 49 |
> |
int nbins = mhisto["Data"]->GetNbinsX(); |
| 50 |
> |
double xmin = mhisto["Data"]->GetXaxis()->GetBinLowEdge(1); |
| 51 |
> |
double xmax = mhisto["Data"]->GetXaxis()->GetBinUpEdge(nbins); |
| 52 |
|
TH1F *test = new TH1F("test", "", nbins, xmin, xmax); |
| 53 |
< |
test -> Sumw2(); |
| 54 |
< |
test -> Add(s1,s2,sf1,sf2); |
| 53 |
> |
test->Sumw2(); |
| 54 |
> |
test->Add(mhisto["QCD"],mhisto["WJets"],sf1,sf2); |
| 55 |
> |
for (map<string,TH1F*>::iterator ihisto = mhisto.begin(); |
| 56 |
> |
ihisto != mhisto.end(); ++ihisto) { |
| 57 |
> |
if (ihisto->first != "QCD" |
| 58 |
> |
&& ihisto->first != "Data" |
| 59 |
> |
&& ihisto->first != "WJets" |
| 60 |
> |
) { |
| 61 |
> |
test->Add(ihisto->second,1.); |
| 62 |
> |
} |
| 63 |
> |
} |
| 64 |
> |
|
| 65 |
|
//test->Scale(1./(1.*test->Integral())); |
| 66 |
< |
Double_t probability = mixS -> KolmogorovTest(test,""); |
| 66 |
> |
Double_t probability = mhisto["Data"]->KolmogorovTest(test,""); |
| 67 |
|
testMC temp = testMC(probability,sf1,sf2); |
| 68 |
|
output.push_back(temp); |
| 69 |
< |
// cout << "probability = " << setw(15) << temp.prob |
| 70 |
< |
// << "; sfQCD = " << setw(10) << temp.scaleF_backg |
| 71 |
< |
// << "; sfWjets = " << setw(6) << temp.scaleF_sample << endl; |
| 69 |
> |
/* cout << "probability = " << setw(15) << temp.prob */ |
| 70 |
> |
/* << "; sfQCD = " << setw(10) << temp.scaleF_backg */ |
| 71 |
> |
/* << "; sfWjets = " << setw(6) << temp.scaleF_sample << endl; */ |
| 72 |
|
delete test; |
| 73 |
|
sf2 = sf2 + stepsize; |
| 74 |
|
} while(sf1 > 0 && sf2 <= 2.0); |
| 76 |
|
} |
| 77 |
|
|
| 78 |
|
//get the maximum KS test result |
| 79 |
< |
testMC getMax(vector<testMC> vec) { |
| 79 |
> |
testMC getMax(vector<testMC> & vec) { |
| 80 |
|
testMC maxKSRes; |
| 81 |
|
Double_t maximum = 0.0; |
| 82 |
+ |
|
| 83 |
+ |
vector<int> posMax; |
| 84 |
|
for (size_t i = 0; i < vec.size(); i++) { |
| 85 |
< |
if (maximum < vec.at(i).prob) { |
| 86 |
< |
maximum = vec.at(i).prob; |
| 87 |
< |
maxKSRes = vec.at(i); |
| 85 |
> |
if (vec.at(i).prob != 1.) { |
| 86 |
> |
if (maximum < vec.at(i).prob) { |
| 87 |
> |
maximum = vec.at(i).prob; |
| 88 |
> |
maxKSRes = vec.at(i); |
| 89 |
> |
} |
| 90 |
> |
} |
| 91 |
> |
else if (vec.at(i).prob == 1.) { |
| 92 |
> |
maximum = 1.; |
| 93 |
> |
posMax.push_back(i); |
| 94 |
> |
} |
| 95 |
> |
} |
| 96 |
> |
|
| 97 |
> |
if (posMax.size() > 0) { |
| 98 |
> |
int tot_ = 0; |
| 99 |
> |
for (vector<int>::iterator iv = posMax.begin(); |
| 100 |
> |
iv != posMax.end(); ++iv) { |
| 101 |
> |
tot_ += *iv; |
| 102 |
|
} |
| 103 |
+ |
int med_ = (int) tot_/posMax.size(); |
| 104 |
+ |
maxKSRes = vec.at(med_); |
| 105 |
|
} |
| 106 |
|
cout << "for maximum: " << setw(12) << maxKSRes.prob |
| 107 |
|
<< "; sb = " << setw(10) << maxKSRes.scaleF_backg |
| 108 |
|
<< "; ss = " << setw(5) << maxKSRes.scaleF_sample << endl; |
| 109 |
|
return maxKSRes; |
| 110 |
|
} |
| 111 |
+ |
|
| 112 |
+ |
Int_t findFirstBinAbove(Double_t threshold,TH1* f) { |
| 113 |
+ |
Int_t nbins = f->GetNbinsX(); |
| 114 |
+ |
for (int nb = 1; nb <= nbins; nb++) { |
| 115 |
+ |
if (f->GetBinContent(nb) > threshold) return nb; |
| 116 |
+ |
} |
| 117 |
+ |
return -1; |
| 118 |
+ |
} |
| 119 |
+ |
|
| 120 |
+ |
Int_t findLastBinAbove(Double_t threshold,TH1* f) { |
| 121 |
+ |
Int_t nbins = f->GetNbinsX(); |
| 122 |
+ |
for (int nb = nbins; nb >= 1; nb--) { |
| 123 |
+ |
if (f->GetBinContent(nb) > threshold) return nb; |
| 124 |
+ |
} |
| 125 |
+ |
return -1; |
| 126 |
+ |
} |
