| 1 |
#include <vector>
|
| 2 |
#include <TFile.h>
|
| 3 |
#include <TH1.h>
|
| 4 |
#include <TH2.h>
|
| 5 |
#include <TLatex.h>
|
| 6 |
#include <TCanvas.h>
|
| 7 |
#include <TMath.h>
|
| 8 |
#include <TLegend.h>
|
| 9 |
#include <TObject.h>
|
| 10 |
#include <iostream>
|
| 11 |
#include <fstream>
|
| 12 |
#include <iomanip>
|
| 13 |
|
| 14 |
using namespace std;
|
| 15 |
|
| 16 |
//define the constants
|
| 17 |
double const weight_QCD = 0.015292368;
|
| 18 |
double const weight_Wjets = 0.002665824;
|
| 19 |
double const weight_ttjets = 0.00009072;
|
| 20 |
Double_t const procQCD = 1.51;
|
| 21 |
Double_t const procWjets = 1.05;
|
| 22 |
Double_t const procttjets = 1.0;
|
| 23 |
|
| 24 |
struct testMC {
|
| 25 |
testMC(Double_t p = 0., Double_t sb = 0., Double_t ss = 0.){prob=p; scaleF_backg = sb; scaleF_sample = ss;}
|
| 26 |
Double_t prob;
|
| 27 |
Double_t scaleF_backg;
|
| 28 |
Double_t scaleF_sample;
|
| 29 |
};
|
| 30 |
|
| 31 |
//function for doing KS test
|
| 32 |
vector<testMC> doKStest(Double_t NmixS, Double_t Ns1, Double_t Ns2, TH1F* mixS, TH1F* s1, TH1F* s2) {
|
| 33 |
vector<testMC> output;
|
| 34 |
//define the scale factors
|
| 35 |
Double_t sf1 = 0.0; // QCD
|
| 36 |
Double_t sf2 = 0.0; // Wjets
|
| 37 |
//KS test
|
| 38 |
do {
|
| 39 |
sf1 = (NmixS - Ns2*sf2)/Ns1;
|
| 40 |
if (sf1 < 0) break;
|
| 41 |
//cout << "..........sf1 = " << sf1 << endl;
|
| 42 |
int nbins = mixS->GetNbinsX();
|
| 43 |
double xmin = mixS->GetXaxis()->GetBinLowEdge(1);
|
| 44 |
double xmax = mixS->GetXaxis()->GetBinUpEdge(nbins);
|
| 45 |
TH1F *test = new TH1F("test", "", nbins, xmin, xmax);
|
| 46 |
test -> Sumw2();
|
| 47 |
test -> Add(s1,s2,sf1,sf2);
|
| 48 |
//test->Scale(1./(1.*test->Integral()));
|
| 49 |
Double_t probability = test -> KolmogorovTest(mixS,"N");
|
| 50 |
testMC temp = testMC(probability,sf1,sf2);
|
| 51 |
output.push_back(temp);
|
| 52 |
// cout << "probability = " << setw(15) << temp.prob
|
| 53 |
// << "; sfQCD = " << setw(10) << temp.scaleF_backg
|
| 54 |
// << "; sfWjets = " << setw(6) << temp.scaleF_sample << endl;
|
| 55 |
delete test;
|
| 56 |
sf2 = sf2 + 0.001;
|
| 57 |
} while(sf1 > 0 && sf2 <= 4.0);
|
| 58 |
return output;
|
| 59 |
}
|
| 60 |
|
| 61 |
//get the maximum KS test result
|
| 62 |
testMC getMax(vector<testMC> vec) {
|
| 63 |
testMC maxKSRes;
|
| 64 |
Double_t maximum = 0.0;
|
| 65 |
for (size_t i = 0; i < vec.size(); i++) {
|
| 66 |
if (maximum < vec.at(i).prob) {
|
| 67 |
maximum = vec.at(i).prob;
|
| 68 |
maxKSRes = vec.at(i);
|
| 69 |
}
|
| 70 |
}
|
| 71 |
cout << "for maximum: " << maxKSRes.prob
|
| 72 |
<< "\tsb = " << maxKSRes.scaleF_backg
|
| 73 |
<< "\tss = " << maxKSRes.scaleF_sample << endl;
|
| 74 |
return maxKSRes;
|
| 75 |
}
|
| 76 |
|
| 77 |
void ikstest() {
|
| 78 |
//Style();
|
| 79 |
TLatex *latex = new TLatex();
|
| 80 |
latex->SetNDC();
|
| 81 |
|
| 82 |
ofstream outprint( "ikstest_results.txt" );
|
| 83 |
|
| 84 |
//open the files with histograms
|
| 85 |
TFile *file[3];
|
| 86 |
file[0] = TFile::Open("InclusiveMu15_Spring10_all.root");
|
| 87 |
file[1] = TFile::Open("WJets_madgraph_Spring10_all.root");
|
| 88 |
file[2] = TFile::Open("TTbarJets_madgraph_Spring10_all.root");
|
| 89 |
|
| 90 |
//define histograms and related parameters
|
| 91 |
string histoN[3] = {"h_mu_pt","h_met0","h_mt0"};
|
| 92 |
string hQCD_new[3] = {"h_mu_pt","h_met0","h_mt0"};
|
| 93 |
string histoLabelX[3] = {"p_{T}^{good Muons}", "E_{T}^{#nu}", "m_{T}^{W}"};
|
| 94 |
Int_t xbins[3] = {10,30,30};
|
| 95 |
Double_t xlow[3] = {0.,0.,0.};
|
| 96 |
Double_t xhigh[3] = {100.,200.,200.};
|
| 97 |
string sample[3] = {"QCD","Wjets","ttjets"};
|
| 98 |
|
| 99 |
TH1F* h_[9];
|
| 100 |
TH1F* mixh_[3];
|
| 101 |
TH1F* hQCD_NEW[3];
|
| 102 |
TH1F* hKSres_[3];
|
| 103 |
TH1F* hKSvalues_[3];
|
| 104 |
|
| 105 |
//load the histograms from the root files
|
| 106 |
for (int i = 0; i < 3; i++) {// 3 variables
|
| 107 |
//cout << "file[" << i << "] : " << endl;
|
| 108 |
string nameNewHisto = "mix_"+histoN[i];
|
| 109 |
string nameNewHistoSFKS = "finalSF_"+histoN[i];
|
| 110 |
string nameNewHistoKSvalues = "KSvalues_"+histoLabelX[i];
|
| 111 |
|
| 112 |
mixh_[i] = new TH1F(nameNewHisto.c_str(),"",xbins[i],xlow[i],xhigh[i]);
|
| 113 |
hKSres_[i] = new TH1F(nameNewHistoSFKS.c_str(),"",xbins[i],xlow[i],xhigh[i]);
|
| 114 |
hKSvalues_[i] = new TH1F(nameNewHistoKSvalues.c_str(),"",300, 0.5, 300.5);
|
| 115 |
//cout << TString("demo/"+hQCD_new[i]) << endl;
|
| 116 |
hQCD_NEW[i] = (TH1F*)file[0]->Get(TString("demo/"+hQCD_new[i])); // corresponds to h_[i]; i=0..2 for now
|
| 117 |
hQCD_NEW[i] -> SetName((hQCD_new[i]).c_str());
|
| 118 |
|
| 119 |
mixh_[i] -> Sumw2();
|
| 120 |
hKSres_[i] -> Sumw2();
|
| 121 |
hKSvalues_[i] -> Sumw2();
|
| 122 |
|
| 123 |
for (int ihisto = 0; ihisto < 3; ihisto++) {
|
| 124 |
//cout << "Variable[" << ihisto << "]" << endl;
|
| 125 |
string histo_name = histoN[ihisto]+sample[i];
|
| 126 |
//cout << TString("demo/"+histoN[ihisto]) << endl;
|
| 127 |
h_[i*3+ihisto] = (TH1F*)file[i]->Get(TString("demo/"+histoN[ihisto]));
|
| 128 |
h_[i*3+ihisto] -> SetName(histo_name.c_str());
|
| 129 |
}
|
| 130 |
}
|
| 131 |
|
| 132 |
//create the mixed samples = "data"
|
| 133 |
for (int i = 0; i < 3; i++) {
|
| 134 |
mixh_[i] -> Add(h_[i],h_[i+3], procQCD,procWjets);
|
| 135 |
//mixh_[i] -> Add(mixh_[i],h_[i+6], 1,procttjets);
|
| 136 |
//cout << "histo_name: " << mixh_[0]->GetNbinsX() << endl;
|
| 137 |
}
|
| 138 |
|
| 139 |
//define the weight corrections for each sample
|
| 140 |
double corr_NevQCD = (h_[2]->GetEntries())*weight_QCD;
|
| 141 |
double corr_NevQCD_NEW = (hQCD_NEW[2]->GetEntries())*weight_QCD;
|
| 142 |
double corr_NevWjets = (h_[5]->GetEntries())*weight_Wjets;
|
| 143 |
double corr_Nevttjets = (h_[8]->GetEntries())*weight_ttjets;
|
| 144 |
double corr_Nevmix = procQCD*corr_NevQCD+procWjets*corr_NevWjets;
|
| 145 |
//double corr_Nevmix = procQCD*corr_NevQCD+procWjets*corr_NevWjets+procttjets*corr_Nevttjets;
|
| 146 |
outprint << "Events mix sample = " << corr_Nevmix << endl;
|
| 147 |
outprint << "Events QCD sample = " << corr_NevQCD << endl;
|
| 148 |
outprint << "Events Wjets sample = " << corr_NevWjets << endl;
|
| 149 |
outprint << "Events Nev revIso QCD sample = " << corr_NevQCD_NEW << endl;
|
| 150 |
|
| 151 |
//define the containers for chosen numbers (coressponding to the max KStest result)
|
| 152 |
testMC maxProb[3];
|
| 153 |
|
| 154 |
//define the scale factors calculated using information obtained from all parameters
|
| 155 |
Double_t SFbackg = 0.0;
|
| 156 |
Double_t sumSFbackg = 0.0;
|
| 157 |
Double_t SFsample = 0.0;
|
| 158 |
Double_t sumSFsample = 0.0;
|
| 159 |
Double_t allKS = 0.0;
|
| 160 |
|
| 161 |
//do the KS test by varying the scale factors
|
| 162 |
for (int i = 0; i < 3; i++) { // 3 variables
|
| 163 |
TH1F *data = (TH1F*)mixh_[i]->Clone();
|
| 164 |
data -> SetName("dataClone");
|
| 165 |
//data -> Scale(1./data->Integral());
|
| 166 |
vector<testMC> resultsKS = doKStest(corr_Nevmix, corr_NevQCD_NEW, corr_NevWjets, data, hQCD_NEW[i], h_[i+3]);
|
| 167 |
testMC tksmax = getMax(resultsKS);
|
| 168 |
maxProb[i] = tksmax;
|
| 169 |
outprint << "\nFor the plot " << histoLabelX[i] << " the results are:"<< endl;
|
| 170 |
outprint << "\tmax Probability = " << maxProb[i].prob << endl;
|
| 171 |
outprint << "\tproc_background = " << maxProb[i].scaleF_backg << endl;
|
| 172 |
outprint << "\tproc_sample = " << maxProb[i].scaleF_sample << endl;
|
| 173 |
|
| 174 |
outprint << "\n\tpercent_B of Data = " << maxProb[i].scaleF_backg*corr_NevQCD_NEW*100/corr_Nevmix << endl;
|
| 175 |
outprint << "\tpercent_S of Data = " << maxProb[i].scaleF_sample*corr_NevWjets*100/corr_Nevmix << endl;
|
| 176 |
outprint << "---------------------------" << endl;
|
| 177 |
|
| 178 |
//create the mixed samples with KS test results
|
| 179 |
sumSFbackg += maxProb[i].prob*maxProb[i].scaleF_backg;
|
| 180 |
sumSFsample += maxProb[i].prob*maxProb[i].scaleF_sample;
|
| 181 |
allKS += maxProb[i].prob;
|
| 182 |
|
| 183 |
//fill a histogram with the results from the KS test for each variable
|
| 184 |
for (int jiter = 0; jiter < resultsKS.size(); jiter++) {
|
| 185 |
//cout << "prob = " << resultsKS.at(jiter).prob << endl;
|
| 186 |
hKSvalues_[i]->SetBinContent(jiter,resultsKS.at(jiter).prob);
|
| 187 |
}
|
| 188 |
delete data;
|
| 189 |
}
|
| 190 |
|
| 191 |
//calculate the final scale factors
|
| 192 |
SFbackg = sumSFbackg/allKS;
|
| 193 |
SFsample = sumSFsample/allKS;
|
| 194 |
outprint << "allKS = " << allKS << "\tbackground = " << SFbackg << "\tsample = " << SFsample << endl;
|
| 195 |
outprint << "\tcombined percent_B of Data = " << SFbackg*corr_NevQCD_NEW*100/corr_Nevmix << endl;
|
| 196 |
outprint << "\tcombined percent_S of Data = " << SFsample*corr_NevWjets*100/corr_Nevmix << endl;
|
| 197 |
outprint << "\n" << endl;
|
| 198 |
outprint << "=================================" << endl;
|
| 199 |
outprint << "\n" << endl;
|
| 200 |
|
| 201 |
}
|
