Jeng/tools/ikstest.h

#include <vector>
#include <TROOT.h>
#include <TFile.h>
#include <TH1.h>
#include <TH2.h>
#include <THStack.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TLegend.h>
#include <TObject.h>
#include <TStyle.h>
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <map>
#include <sys/types.h>
#include <sys/stat.h>
#include <stdio.h>

const double stepsize = 0.001;
//=================================
// Structure/Functions
//=================================

struct testMC {
  testMC(Double_t p = 0., Double_t sb = 0., Double_t ss = 0.){prob=p; scaleF_backg = sb; scaleF_sample = ss;}
  Double_t prob;
  Double_t scaleF_backg;
  Double_t scaleF_sample;
};

//function for doing KS test
vector<testMC> doKStest(const vector<Double_t>& vnevt, map<string,TH1F*>& mhisto) {
  vector<testMC> output;
  //define the scale factors
  Double_t sf1 = 0.0; // QCD
  Double_t sf2 = 0.0; // Wjets
  Double_t Ns  = 0.0; // non-Wjets (ttbar,Zjets,single top...)
  for (UInt_t i = 3; i < vnevt.size(); ++i)
    Ns += vnevt[i];
  cout << "Total number of Non Wjets events = " << Ns << endl;
  //KS test
  do {
    sf1 = (vnevt.at(0) - Ns - vnevt.at(2)*sf2)/vnevt.at(1);
    if (sf1 < 0) break;
    //cout << "..........sf1 = " << sf1 << endl;
    int nbins = mhisto["Data"]->GetNbinsX();
    double xmin = mhisto["Data"]->GetXaxis()->GetBinLowEdge(1);
    double xmax = mhisto["Data"]->GetXaxis()->GetBinUpEdge(nbins);
    TH1F *test = new TH1F("test", "", nbins, xmin, xmax);
    test->Sumw2();
    test->Add(mhisto["QCD"],mhisto["WJets"],sf1,sf2);
    for (map<string,TH1F*>::iterator ihisto = mhisto.begin();
         ihisto != mhisto.end(); ++ihisto) {
      if (ihisto->first != "QCD" 
          && ihisto->first != "Data"
          && ihisto->first != "WJets"
          ) {
        test->Add(ihisto->second,1.);
      }
    }

    //test->Scale(1./(1.*test->Integral()));
    Double_t probability = mhisto["Data"]->KolmogorovTest(test,"");
    testMC temp = testMC(probability,sf1,sf2);
    output.push_back(temp);
/*     cout << "probability = " << setw(15) << temp.prob  */
/*       << "; sfQCD = "     << setw(10) << temp.scaleF_backg  */
/*       << "; sfWjets = "   << setw(6)  << temp.scaleF_sample << endl; */
    delete test;
    sf2 = sf2 + stepsize;
  } while(sf1 > 0 && sf2 <= 2.0);
  return output;
}

//get the maximum KS test result
testMC getMax(const vector<testMC> & vec) {
  testMC maxKSRes;
  Double_t maximum = 0.0;

  vector<int> posMax;
  for (size_t i = 0; i < vec.size(); i++) {
    if (vec.at(i).prob != 1.) {
      if (maximum < vec.at(i).prob)  {
        maximum = vec.at(i).prob;
        maxKSRes = vec.at(i);
      }
    }
    else if (vec.at(i).prob == 1.) {
      maximum = 1.;
      posMax.push_back(i);
    }
  }

  if (posMax.size() > 0) {
    int tot_ = 0;
    for (vector<int>::iterator iv = posMax.begin();
         iv != posMax.end(); ++iv) {
      tot_ += *iv;
    }
    int med_ = (int) tot_/posMax.size();
    maxKSRes = vec.at(med_);
  }
  cout << "[KSmax] for maximum: " << setw(12) << maxKSRes.prob 
       << "; sb = " << setw(10) << maxKSRes.scaleF_backg 
       << "; ss = " << setw(5) << maxKSRes.scaleF_sample << endl; 
  return maxKSRes;
}

testMC getMaxL(const vector<testMC> vec[]) {
  testMC maxKSRes;
  Double_t maximum = 0.0;
  vector<int> posMax;

  for (size_t i = 0; i < vec[0].size(); i++) {
    Double_t product_ = 1.0;
    for (size_t j = 0; j < 3; j++) // 3 variables
      product_ *= vec[j].at(i).prob;
    if (product_ != 1.) {
      if (product_ > maximum) {
        maximum = product_;
        maxKSRes = vec[0].at(i);
      }
    }
    else {
      maximum = 1.;
      posMax.push_back(i);
    }
  }

  if (posMax.size() > 0) {
    int tot_ = 0;
    for (vector<int>::iterator iv = posMax.begin();
         iv != posMax.end(); ++iv) {
      tot_ += *iv;
    }
    int med_ = (int) tot_/posMax.size();
    maxKSRes = vec[0].at(med_);
  }
  maxKSRes.prob = maximum;
  cout << "[KSopt] for maximum: " << setw(12) << maxKSRes.prob 
       << "; sb = " << setw(10) << maxKSRes.scaleF_backg 
       << "; ss = " << setw(5) << maxKSRes.scaleF_sample << endl; 
  return maxKSRes;
}

Int_t findFirstBinAbove(Double_t threshold,TH1* f) {
  Int_t nbins = f->GetNbinsX();
  for (int nb = 1; nb <= nbins; nb++) {
    if (f->GetBinContent(nb) > threshold) return nb;
  }
  return -1;
}

Int_t findLastBinAbove(Double_t threshold,TH1* f) {
  Int_t nbins = f->GetNbinsX();
  for (int nb = nbins; nb >= 1; nb--) {
    if (f->GetBinContent(nb) > threshold) return nb;
  }
  return -1;
}
Revision:	1.2
Committed:	Thu Sep 16 04:09:04 2010 UTC (14 years, 7 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	gregj-20110224, gregj-20110223-v1, gregj-20110223, gregj-20110205, gregj-20110114, gregj-20110107, gregj-20101223, gregj-20101212, gregj-20101211, gregj-20101210-new, gregj-20101210, gregj-20101121, gregj-20101113-NEW, gregj-20101113, gregj-20101111-JES_JETONLY, gregj-20101103, gregj-20101101, gregj-20101019, gregj-20101005, greg-20100929, gregj-20100925, gregj-20100916, HEAD
Changes since 1.1:	+39 -2 lines
Log Message:	update
#	User	Rev	Content
1	jengbou	1.1	#include <vector>
2			#include <TROOT.h>
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>
16			#include <iomanip>
17			#include <map>
18			#include <sys/types.h>
19			#include <sys/stat.h>
20			#include <stdio.h>
21
22			const double stepsize = 0.001;
23			//=================================
24			// Structure/Functions
25			//=================================
26
27			struct testMC {
28			testMC(Double_t p = 0., Double_t sb = 0., Double_t ss = 0.){prob=p; scaleF_backg = sb; scaleF_sample = ss;}
29			Double_t prob;
30			Double_t scaleF_backg;
31			Double_t scaleF_sample;
32			};
33
34			//function for doing KS test
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 = (vnevt.at(0) - Ns - vnevt.at(2)*sf2)/vnevt.at(1);
47			if (sf1 < 0) break;
48			//cout << "..........sf1 = " << sf1 << endl;
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(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 = 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; */
72			delete test;
73			sf2 = sf2 + stepsize;
74			} while(sf1 > 0 && sf2 <= 2.0);
75			return output;
76			}
77
78			//get the maximum KS test result
79	jengbou	1.2	testMC getMax(const vector<testMC> & vec) {
80	jengbou	1.1	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 (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	jengbou	1.2	cout << "[KSmax] 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			testMC getMaxL(const vector<testMC> vec[]) {
113			testMC maxKSRes;
114			Double_t maximum = 0.0;
115			vector<int> posMax;
116
117			for (size_t i = 0; i < vec[0].size(); i++) {
118			Double_t product_ = 1.0;
119			for (size_t j = 0; j < 3; j++) // 3 variables
120			product_ *= vec[j].at(i).prob;
121			if (product_ != 1.) {
122			if (product_ > maximum) {
123			maximum = product_;
124			maxKSRes = vec[0].at(i);
125			}
126			}
127			else {
128			maximum = 1.;
129			posMax.push_back(i);
130			}
131			}
132
133			if (posMax.size() > 0) {
134			int tot_ = 0;
135			for (vector<int>::iterator iv = posMax.begin();
136			iv != posMax.end(); ++iv) {
137			tot_ += *iv;
138			}
139			int med_ = (int) tot_/posMax.size();
140			maxKSRes = vec[0].at(med_);
141			}
142			maxKSRes.prob = maximum;
143			cout << "[KSopt] for maximum: " << setw(12) << maxKSRes.prob
144	jengbou	1.1	<< "; sb = " << setw(10) << maxKSRes.scaleF_backg
145			<< "; ss = " << setw(5) << maxKSRes.scaleF_sample << endl;
146			return maxKSRes;
147			}
148
149			Int_t findFirstBinAbove(Double_t threshold,TH1* f) {
150			Int_t nbins = f->GetNbinsX();
151			for (int nb = 1; nb <= nbins; nb++) {
152			if (f->GetBinContent(nb) > threshold) return nb;
153			}
154			return -1;
155			}
156
157			Int_t findLastBinAbove(Double_t threshold,TH1* f) {
158			Int_t nbins = f->GetNbinsX();
159			for (int nb = nbins; nb >= 1; nb--) {
160			if (f->GetBinContent(nb) > threshold) return nb;
161			}
162			return -1;
163			}