macros/src/MatrixMethod_jetbin_separately.cc

#include "TFile.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH2.h"
#include "TH2D.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TLatex.h"

#include "../../../Style.cc"
#include "MatrixMethod.h"

using namespace std;

void mainMM() {
  Style();
  TLatex *latex = new TLatex();
  latex->SetNDC();
    
  string path = "/uscms/home/omanghel/work/Top/macros/";
  string dirName = "FastSimWinter09";
  string fileName = "histosforMM_revd0cut_2jet.root";
  string txtfileName = "inputMM_true_2jet.txt";
  //string cvtxtfileName = "inputMM_.txt";
  
  //define the variable to store the statistical information
  string sample[5];
  double weightf[5];
  double NlooseTrue[5];
  double NtightTrue[5];
  double NcontrolVarloose[5];
  double NcontrolVartight[5];
  int Nsamples = 5;
  for (int i=0; i< Nsamples; i++) {
        weightf[i] = 0.0;
        NlooseTrue[i] = 0;
        NtightTrue[i] = 0;
        NcontrolVarloose[i] = 0;
        NcontrolVartight[i] = 0;
  }
        
  //read the statistics from the ascii file
  ifstream txtFile;
  txtFile.open(txtfileName.c_str());
  string temp0; 
  double temp1;
  int temp2;
  int temp3;
  int temp4;
  int temp5;
  for (int i=0; i < Nsamples; i++) {
          txtFile >> temp0 >> temp1 >> temp2 >> temp3 >> temp4 >> temp5;
          sample[i] = temp0;
          weightf[i] = temp1;
          NlooseTrue[i] = temp2;
          NtightTrue[i] = temp3;
          NcontrolVarloose[i] = temp4;
          NcontrolVartight[i] = temp5;
  }
  txtFile.close();
        
  //add the weighted events to form a data sample
  double Nloose_data = 0.0;
  double Ntight_data = 0.0;
  double Nsigloose_MC = 0.0;
  double Nsigtight_MC = 0.0;
  double NcontrolVarDataloose = 0.0;
  double NcontrolVarDatatight = 0.0;
  for (int i=0; i != Nsamples; i++) {
          cout << sample[i] << "\t" << weightf[i] << "\t" << NlooseTrue[i] << "\t" << NtightTrue[i] << "\t" << NcontrolVarloose[i] << "\t" << NcontrolVartight[i] << endl;
        Nloose_data += NlooseTrue[i]*weightf[i];
        Ntight_data += NtightTrue[i]*weightf[i];
        NcontrolVarDataloose += NcontrolVarloose[i]*weightf[i];
        NcontrolVarDatatight += NcontrolVartight[i]*weightf[i];
        if (i!=0) {
                Nsigloose_MC += NlooseTrue[i]*weightf[i];
                Nsigtight_MC += NtightTrue[i]*weightf[i];
        }
  }
  double NLminT_data = Nloose_data - Ntight_data;
        
  //get the corresponding tagging efficiencies and their errors
  double epsSig = Nsigtight_MC*1.0/Nsigloose_MC;
  double DepsSig = 0.05; //based only on weight factors approximations
  //double epsQCD = NcontrolVartight[0]*1.0/NcontrolVarloose[0];
  double epsQCD = NcontrolVarDatatight*1.0/NcontrolVarDataloose;
  double DepsQCD = 0.002 ; //based on small variations around the cut value
        
  cout << "epsQCD = " << epsQCD << endl;
  cout << "epsSig = " << epsSig << endl;

  //apply the matrix method
  outInfo result = mxCalculate(NLminT_data,Ntight_data,epsSig,DepsSig,epsQCD,DepsQCD);
        
  cout << "NSig = " << result.nSig << endl;
  cout << "ErrSig = " << result.errSig << endl;
  cout << "NBkg = " << result.nBkg << endl;
  cout << "ErrBkg = " << result.errBkg << endl; 
  cout << "NSig_true = " << Nsigtight_MC << endl;
  cout << "NBkg_true = " << NtightTrue[0]*weightf[0] << endl;

 ///////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////
 /////////////////Plot the histograms//////////////////////////
 //////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////
        
  TFile *inFile;
  //inFile = TFile::Open((path+"/"+dirName+"/"+fileName).c_str());
  inFile = TFile::Open(fileName.c_str());
  inFile->cd();
        
  //string sample[5] = {"inclusiveMu_none","TTjets_none","Wjets_Wbb","Wjets_Wcc","Wjets_Wjets"};
  string selection = "loose";
  string histo1DType[6] = {"d0d0sig","dPhi_MuandMet","dPhi_Muonand2Jets","dPhi_MuonandJet","metEt","muPt"};
  string xlabel[6] = {"d0/d0sig","|#Delta#phi(#mu,MET)|","|#Delta#phi(#mu,jet1)|+|#Delta#phi(#mu,jet2)|","|#Delta#phi(#mu,jet)|","E_{T}^{met}","p_{T}^{#mu}"};
  string ylabel[6] = {"Entries","Entries","Entries","Entries","Entries","Entries"};
  string histo2DType[4] = {"muRelIso_vs_d0corr","muRelIso_vs_dPhimuMet","muRelIso_vs_metEt","muRelIso_vs_muPt"};
  string xlabel2D[4] = {"Isolation^{#mu}","Isolation^{#mu}","Isolation^{#mu}","Isolation^{#mu}"};
  string ylabel2D[4] = {"d0/d0sig","|#Delta#phi(#mu,MET)|","E_{T}^{met}","p_{T}^{#mu}"};
        
  Double_t minBin[6] = {0,0,0,0,0,0};
  Double_t maxBin[6] = {3,4,10,4,500,500};
  int color[5] = {1,2,3,4,6};  //for the Nsamples
  
  TH1F* histo1D[30];
  TH2D* histo2D[20]; //for the Nsamples
  
  for (int i=0; i!=6; i++) {
        for (int j=0; j != Nsamples; j++) {
                string nametmp;
                if (j < 2) nametmp = "h_"+histo1DType[i]+"_"+selection+"_"+sample[j]+"nod0cut";
                else nametmp = "h_"+histo1DType[i]+"_"+selection+"_"+sample[j];
                cout << nametmp << endl;
                gDirectory->GetObject(nametmp.c_str(),histo1D[i*Nsamples+j]);
                //histo1D[i*Nsamples+j] = inFile->Get(nametmp.c_str());
                histo1D[i*Nsamples+j]->Scale(1/histo1D[i*Nsamples+j]->Integral());
                histo1D[i*Nsamples+j]->SetLineColor(color[j]);
                histo1D[i*Nsamples+j]->SetLineWidth(2.0);
                histo1D[i*Nsamples+j]->GetXaxis()->SetTitle(xlabel[i].c_str());
                histo1D[i*Nsamples+j]->GetYaxis()->SetTitle(ylabel[i].c_str());
                if (i==0) histo1D[i*Nsamples+j]->GetYaxis()->SetRangeUser(0.0,0.09);
                if (i==2) histo1D[i*Nsamples+j]->GetYaxis()->SetRangeUser(0.0,0.055);
                histo1D[i*Nsamples+j]->SetStats(0);
        
                if (i < 4) {
                        string nametmp2D; 
                        if (j < 2) nametmp2D = "h_"+histo2DType[i]+"_"+selection+"_"+sample[j]+"nod0cut";
                        else nametmp2D = "h_"+histo2DType[i]+"_"+selection+"_"+sample[j];
                        cout << nametmp2D << endl;
                        //if (i==0) {
                                gDirectory->GetObject(nametmp2D.c_str(),histo2D[i*Nsamples+j]);
                                histo2D[i*Nsamples+j]->SetLineColor(color[j]);
                                histo2D[i*Nsamples+j]->SetMarkerColor(color[j]);
                                histo2D[i*Nsamples+j]->GetXaxis()->SetTitle(xlabel2D[i].c_str());
                                histo2D[i*Nsamples+j]->GetYaxis()->SetTitle(ylabel2D[i].c_str());
                        //}
                }
        }
        
        TCanvas *canvas1 = new TCanvas(histo1DType[i].c_str(), "");
                cout << "here4" << endl;
        histo1D[i*Nsamples]->Draw("");
                cout << "here5" << endl;
        for (int j=0; j != Nsamples; j++) {
                if (j==0) continue;
                histo1D[i*Nsamples+j]->Draw("same");
                cout << "here6" << endl;
                }
        TLegend *legend1 = new TLegend(0.6, 0.6, 0.9, 0.85);
        for (int j=0; j < Nsamples; j++) legend1->AddEntry(histo1D[i*Nsamples+j],sample[j].c_str());
        legend1->Draw();
        legend1->SetFillColor(kWhite);
        latex->DrawLatex(0.22,0.91,histo1DType[i].c_str());
        //canvas1->SetLogy();
        string nameCanvas1 = histo1DType[i]+".eps";
        canvas1->SaveAs(nameCanvas1.c_str());

                cout << "here7" << endl;
        if (i < 4) {
                //new TCanvas();
                TCanvas *canvas2 = new TCanvas(histo2DType[i].c_str(), "");
                histo2D[i*Nsamples]->Draw("");
                for (int k=0; k!=Nsamples; k++) {
                                if (k==0) continue;
                        histo2D[i*Nsamples+k]->Draw("same");
                        }
                TLegend *legend2 = new TLegend(0.7, 0.5, 0.9, 0.85);
                for (int k=0; k < Nsamples; k++) legend2->AddEntry(histo2D[i*Nsamples+k],sample[k].c_str());
                legend2->Draw();
                legend2->SetFillColor(kWhite);
                latex->DrawLatex(0.22,0.91,histo2DType[i].c_str());
                //canvas1->SetLogy();
                string nameCanvas2 = histo2DType[i]+".eps";
                canvas2->SaveAs(nameCanvas2.c_str());
                //cout << "here8"<< endl;
        }       
  }
}
Revision:	1.1
Committed:	Fri Jan 29 23:04:55 2010 UTC (15 years, 3 months ago) by anghel
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	Revision 29Jan2010
#	User	Rev	Content
1	anghel	1.1	#include "TFile.h"
2			#include "TH1.h"
3			#include "TH1F.h"
4			#include "TH2.h"
5			#include "TH2D.h"
6			#include "TFile.h"
7			#include "TDirectory.h"
8			#include "TCanvas.h"
9			#include "TLegend.h"
10			#include "TLatex.h"
11
12			#include "../../../Style.cc"
13			#include "MatrixMethod.h"
14
15			using namespace std;
16
17			void mainMM() {
18			Style();
19			TLatex *latex = new TLatex();
20			latex->SetNDC();
21
22			string path = "/uscms/home/omanghel/work/Top/macros/";
23			string dirName = "FastSimWinter09";
24			string fileName = "histosforMM_revd0cut_2jet.root";
25			string txtfileName = "inputMM_true_2jet.txt";
26			//string cvtxtfileName = "inputMM_.txt";
27
28			//define the variable to store the statistical information
29			string sample[5];
30			double weightf[5];
31			double NlooseTrue[5];
32			double NtightTrue[5];
33			double NcontrolVarloose[5];
34			double NcontrolVartight[5];
35			int Nsamples = 5;
36			for (int i=0; i< Nsamples; i++) {
37			weightf[i] = 0.0;
38			NlooseTrue[i] = 0;
39			NtightTrue[i] = 0;
40			NcontrolVarloose[i] = 0;
41			NcontrolVartight[i] = 0;
42			}
43
44			//read the statistics from the ascii file
45			ifstream txtFile;
46			txtFile.open(txtfileName.c_str());
47			string temp0;
48			double temp1;
49			int temp2;
50			int temp3;
51			int temp4;
52			int temp5;
53			for (int i=0; i < Nsamples; i++) {
54			txtFile >> temp0 >> temp1 >> temp2 >> temp3 >> temp4 >> temp5;
55			sample[i] = temp0;
56			weightf[i] = temp1;
57			NlooseTrue[i] = temp2;
58			NtightTrue[i] = temp3;
59			NcontrolVarloose[i] = temp4;
60			NcontrolVartight[i] = temp5;
61			}
62			txtFile.close();
63
64			//add the weighted events to form a data sample
65			double Nloose_data = 0.0;
66			double Ntight_data = 0.0;
67			double Nsigloose_MC = 0.0;
68			double Nsigtight_MC = 0.0;
69			double NcontrolVarDataloose = 0.0;
70			double NcontrolVarDatatight = 0.0;
71			for (int i=0; i != Nsamples; i++) {
72			cout << sample[i] << "\t" << weightf[i] << "\t" << NlooseTrue[i] << "\t" << NtightTrue[i] << "\t" << NcontrolVarloose[i] << "\t" << NcontrolVartight[i] << endl;
73			Nloose_data += NlooseTrue[i]*weightf[i];
74			Ntight_data += NtightTrue[i]*weightf[i];
75			NcontrolVarDataloose += NcontrolVarloose[i]*weightf[i];
76			NcontrolVarDatatight += NcontrolVartight[i]*weightf[i];
77			if (i!=0) {
78			Nsigloose_MC += NlooseTrue[i]*weightf[i];
79			Nsigtight_MC += NtightTrue[i]*weightf[i];
80			}
81			}
82			double NLminT_data = Nloose_data - Ntight_data;
83
84			//get the corresponding tagging efficiencies and their errors
85			double epsSig = Nsigtight_MC*1.0/Nsigloose_MC;
86			double DepsSig = 0.05; //based only on weight factors approximations
87			//double epsQCD = NcontrolVartight[0]*1.0/NcontrolVarloose[0];
88			double epsQCD = NcontrolVarDatatight*1.0/NcontrolVarDataloose;
89			double DepsQCD = 0.002 ; //based on small variations around the cut value
90
91			cout << "epsQCD = " << epsQCD << endl;
92			cout << "epsSig = " << epsSig << endl;
93
94			//apply the matrix method
95			outInfo result = mxCalculate(NLminT_data,Ntight_data,epsSig,DepsSig,epsQCD,DepsQCD);
96
97			cout << "NSig = " << result.nSig << endl;
98			cout << "ErrSig = " << result.errSig << endl;
99			cout << "NBkg = " << result.nBkg << endl;
100			cout << "ErrBkg = " << result.errBkg << endl;
101			cout << "NSig_true = " << Nsigtight_MC << endl;
102			cout << "NBkg_true = " << NtightTrue[0]*weightf[0] << endl;
103
104			///////////////////////////////////////////////////////////////
105			///////////////////////////////////////////////////////////////
106			/////////////////Plot the histograms//////////////////////////
107			//////////////////////////////////////////////////////////////
108			///////////////////////////////////////////////////////////////
109
110			TFile *inFile;
111			//inFile = TFile::Open((path+"/"+dirName+"/"+fileName).c_str());
112			inFile = TFile::Open(fileName.c_str());
113			inFile->cd();
114
115			//string sample[5] = {"inclusiveMu_none","TTjets_none","Wjets_Wbb","Wjets_Wcc","Wjets_Wjets"};
116			string selection = "loose";
117			string histo1DType[6] = {"d0d0sig","dPhi_MuandMet","dPhi_Muonand2Jets","dPhi_MuonandJet","metEt","muPt"};
118			string xlabel[6] = {"d0/d0sig","\|#Delta#phi(#mu,MET)\|","\|#Delta#phi(#mu,jet1)\|+\|#Delta#phi(#mu,jet2)\|","\|#Delta#phi(#mu,jet)\|","E_{T}^{met}","p_{T}^{#mu}"};
119			string ylabel[6] = {"Entries","Entries","Entries","Entries","Entries","Entries"};
120			string histo2DType[4] = {"muRelIso_vs_d0corr","muRelIso_vs_dPhimuMet","muRelIso_vs_metEt","muRelIso_vs_muPt"};
121			string xlabel2D[4] = {"Isolation^{#mu}","Isolation^{#mu}","Isolation^{#mu}","Isolation^{#mu}"};
122			string ylabel2D[4] = {"d0/d0sig","\|#Delta#phi(#mu,MET)\|","E_{T}^{met}","p_{T}^{#mu}"};
123
124			Double_t minBin[6] = {0,0,0,0,0,0};
125			Double_t maxBin[6] = {3,4,10,4,500,500};
126			int color[5] = {1,2,3,4,6}; //for the Nsamples
127
128			TH1F* histo1D[30];
129			TH2D* histo2D[20]; //for the Nsamples
130
131			for (int i=0; i!=6; i++) {
132			for (int j=0; j != Nsamples; j++) {
133			string nametmp;
134			if (j < 2) nametmp = "h_"+histo1DType[i]+"_"+selection+"_"+sample[j]+"nod0cut";
135			else nametmp = "h_"+histo1DType[i]+"_"+selection+"_"+sample[j];
136			cout << nametmp << endl;
137			gDirectory->GetObject(nametmp.c_str(),histo1D[i*Nsamples+j]);
138			//histo1D[i*Nsamples+j] = inFile->Get(nametmp.c_str());
139			histo1D[iNsamples+j]->Scale(1/histo1D[iNsamples+j]->Integral());
140			histo1D[i*Nsamples+j]->SetLineColor(color[j]);
141			histo1D[i*Nsamples+j]->SetLineWidth(2.0);
142			histo1D[i*Nsamples+j]->GetXaxis()->SetTitle(xlabel[i].c_str());
143			histo1D[i*Nsamples+j]->GetYaxis()->SetTitle(ylabel[i].c_str());
144			if (i==0) histo1D[i*Nsamples+j]->GetYaxis()->SetRangeUser(0.0,0.09);
145			if (i==2) histo1D[i*Nsamples+j]->GetYaxis()->SetRangeUser(0.0,0.055);
146			histo1D[i*Nsamples+j]->SetStats(0);
147
148			if (i < 4) {
149			string nametmp2D;
150			if (j < 2) nametmp2D = "h_"+histo2DType[i]+"_"+selection+"_"+sample[j]+"nod0cut";
151			else nametmp2D = "h_"+histo2DType[i]+"_"+selection+"_"+sample[j];
152			cout << nametmp2D << endl;
153			//if (i==0) {
154			gDirectory->GetObject(nametmp2D.c_str(),histo2D[i*Nsamples+j]);
155			histo2D[i*Nsamples+j]->SetLineColor(color[j]);
156			histo2D[i*Nsamples+j]->SetMarkerColor(color[j]);
157			histo2D[i*Nsamples+j]->GetXaxis()->SetTitle(xlabel2D[i].c_str());
158			histo2D[i*Nsamples+j]->GetYaxis()->SetTitle(ylabel2D[i].c_str());
159			//}
160			}
161			}
162
163			TCanvas *canvas1 = new TCanvas(histo1DType[i].c_str(), "");
164			cout << "here4" << endl;
165			histo1D[i*Nsamples]->Draw("");
166			cout << "here5" << endl;
167			for (int j=0; j != Nsamples; j++) {
168			if (j==0) continue;
169			histo1D[i*Nsamples+j]->Draw("same");
170			cout << "here6" << endl;
171			}
172			TLegend *legend1 = new TLegend(0.6, 0.6, 0.9, 0.85);
173			for (int j=0; j < Nsamples; j++) legend1->AddEntry(histo1D[i*Nsamples+j],sample[j].c_str());
174			legend1->Draw();
175			legend1->SetFillColor(kWhite);
176			latex->DrawLatex(0.22,0.91,histo1DType[i].c_str());
177			//canvas1->SetLogy();
178			string nameCanvas1 = histo1DType[i]+".eps";
179			canvas1->SaveAs(nameCanvas1.c_str());
180
181			cout << "here7" << endl;
182			if (i < 4) {
183			//new TCanvas();
184			TCanvas *canvas2 = new TCanvas(histo2DType[i].c_str(), "");
185			histo2D[i*Nsamples]->Draw("");
186			for (int k=0; k!=Nsamples; k++) {
187			if (k==0) continue;
188			histo2D[i*Nsamples+k]->Draw("same");
189			}
190			TLegend *legend2 = new TLegend(0.7, 0.5, 0.9, 0.85);
191			for (int k=0; k < Nsamples; k++) legend2->AddEntry(histo2D[i*Nsamples+k],sample[k].c_str());
192			legend2->Draw();
193			legend2->SetFillColor(kWhite);
194			latex->DrawLatex(0.22,0.91,histo2DType[i].c_str());
195			//canvas1->SetLogy();
196			string nameCanvas2 = histo2DType[i]+".eps";
197			canvas2->SaveAs(nameCanvas2.c_str());
198			//cout << "here8"<< endl;
199			}
200			}
201			}