Jeng/scripts/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(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 << "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.4
Committed:	Thu Sep 9 08:03:25 2010 UTC (14 years, 8 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.3:	+40 -5 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	jengbou	1.2	#include <THStack.h>
7	jengbou	1.1	#include <TLatex.h>
8			#include <TCanvas.h>
9			#include <TMath.h>
10			#include <TLegend.h>
11			#include <TObject.h>
12	jengbou	1.4	#include <TStyle.h>
13	jengbou	1.1	#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	jengbou	1.3	vector<testMC> doKStest(const vector<Double_t>& vnevt, map<string,TH1F*>& mhisto) {
36	jengbou	1.1	vector<testMC> output;
37			//define the scale factors
38			Double_t sf1 = 0.0; // QCD
39			Double_t sf2 = 0.0; // Wjets
40	jengbou	1.3	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	jengbou	1.4	cout << "Total number of Non Wjets events = " << Ns << endl;
44	jengbou	1.1	//KS test
45			do {
46	jengbou	1.3	sf1 = (vnevt.at(0) - Ns - vnevt.at(2)*sf2)/vnevt.at(1);
47	jengbou	1.1	if (sf1 < 0) break;
48			//cout << "..........sf1 = " << sf1 << endl;
49	jengbou	1.3	int nbins = mhisto["Data"]->GetNbinsX();
50			double xmin = mhisto["Data"]->GetXaxis()->GetBinLowEdge(1);
51			double xmax = mhisto["Data"]->GetXaxis()->GetBinUpEdge(nbins);
52	jengbou	1.1	TH1F *test = new TH1F("test", "", nbins, xmin, xmax);
53	jengbou	1.3	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	jengbou	1.1	//test->Scale(1./(1.*test->Integral()));
66	jengbou	1.3	Double_t probability = mhisto["Data"]->KolmogorovTest(test,"");
67	jengbou	1.1	testMC temp = testMC(probability,sf1,sf2);
68			output.push_back(temp);
69	jengbou	1.3	/* cout << "probability = " << setw(15) << temp.prob */
70			/* << "; sfQCD = " << setw(10) << temp.scaleF_backg */
71			/* << "; sfWjets = " << setw(6) << temp.scaleF_sample << endl; */
72	jengbou	1.1	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.4	testMC getMax(vector<testMC> & vec) {
80	jengbou	1.1	testMC maxKSRes;
81			Double_t maximum = 0.0;
82	jengbou	1.4
83			vector<int> posMax;
84	jengbou	1.1	for (size_t i = 0; i < vec.size(); i++) {
85	jengbou	1.4	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	jengbou	1.1	}
103	jengbou	1.4	int med_ = (int) tot_/posMax.size();
104			maxKSRes = vec.at(med_);
105	jengbou	1.1	}
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	jengbou	1.4
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			}