MitHzz4l/Shapes/makeShape.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TROOT.h>                   // access to gROOT, entry point to ROOT system
#include <TLegend.h>
#include <TSystem.h>                 // interface to OS
#include <TFile.h>                   // file handle class
#include <TCanvas.h>
#include <TH1F.h>
#include <TTree.h>                   // class to access ntuples
#include <TClonesArray.h>            // ROOT array class
#include <TLorentzVector.h>           // 4-vector class
#include <TBenchmark.h>              // class to track macro running statistics
#include <iostream>                  // standard I/O
#include <iomanip>                   // functions to format standard I/O
#include <fstream>                   // functions for file I/O
#include <sstream>                   // class for parsing strings
#include <map>

#endif

TCanvas *c1 = new TCanvas("c1","c1");
vector<float> masses, widths;
map <TString,TString> binstrings;
//----------------------------------------------------------------------------------------
void formatBins()
{
  binstrings["def"] = "(75,30,500)";

  masses.push_back(115);
  masses.push_back(120);
  masses.push_back(130);
  masses.push_back(140);
  masses.push_back(150);
  masses.push_back(160);
  masses.push_back(170);  
  masses.push_back(180);
  masses.push_back(190);
  masses.push_back(200);
  masses.push_back(230);
  masses.push_back(250);
  masses.push_back(300);
  masses.push_back(350);
  masses.push_back(400);
  masses.push_back(500);
  masses.push_back(600);

  widths.push_back(0.0166566);
  widths.push_back(0.0185038);
  widths.push_back(0.0259764);
  widths.push_back(0.0425892);
  widths.push_back(0.0920884);
  widths.push_back(0.438603 );
  widths.push_back(1.96908  );
  widths.push_back(3.30333  );
  widths.push_back(5.40595  );
  widths.push_back(7.42329  );
  widths.push_back(14.688   );
  widths.push_back(20.9409  );
  widths.push_back(43.5433  );
  widths.push_back(77.7866  );
  widths.push_back(141.866  );
  widths.push_back(266.492  );
  widths.push_back(379.294  );

  stringstream ss;
  for(unsigned i=0; i<masses.size(); i++) {
    float mass = masses[i];
    float width = widths[i];
    ss.str("");
    ss << "(30," << mass-width << "," << mass+width << ")";
    TString binning(ss.str());
    ss.str("");
    ss << mass;
    binstrings[ss.str()] = binning;
  }
}
//----------------------------------------------------------------------------------------
TH1F *getMCFM(TString mass, TString process, TString production, TString updown)
{
  TString scale;
  if(updown=="def")     scale = "91_";
  else if(updown=="lo") scale = "45_";
  else if(updown=="hi") scale = "182";
  else assert(0);

  TString fname,wgtStr;
  if(process=="higgs")   {
    fname = "/temp/dkralph/mcfm/114-"+mass+"-"+updown+"/HZZ_4l_tota_mstw8nl_"+scale+"_"+scale+"_.root";
    wgtStr = "wt_gg";
  } else if(process=="zz") {
    fname = "/temp/dkralph/mcfm/81-115-"+updown+"/ZZlept_tota_mstw8nl_"+scale+"_"+scale+"_.root";
    if(production=="gg")          wgtStr = "wt_gg";
    else if(production=="other")  wgtStr = "wt_gq+wt_qqb";
    else assert(0);
  } else assert(0);

  TFile f(fname); assert(f.IsOpen());
  TTree *tree=0;
  f.GetObject("h10",tree); assert(tree);
  TString A,B,C,D,m,varexp;
  A = "(px3+px4+px5+px6)*(px3+px4+px5+px6)";
  B = "(py3+py4+py5+py6)*(py3+py4+py5+py6)";
  C = "(pz3+pz4+pz5+pz6)*(pz3+pz4+pz5+pz6)";
  D = "(E_3+E_4+E_5+E_6)*(E_3+E_4+E_5+E_6)";
  m = "sqrt("+D+"-"+A+"-"+B+"-"+C+")";
  TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
  varexp= m+">>m4l"+binning;
  tree->Draw(varexp,wgtStr,"",500000);
  TH1F *hm = (TH1F*)gDirectory->Get("m4l");
  hm->SetDirectory(0);
  hm->SetTitle("m4l");
  f.Close();
  hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
  return hm;
}  
//----------------------------------------------------------------------------------------
TH1F *getPowheg(TString mass, TString process, TString production)
{
  TString fname;
  if(process=="higgs") {
    if(production=="gg") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu.root";
    else assert(0);
  } else if(process=="zz") {
    if(production=="gg")         fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-gg2zz_ntuple.root";
    else if(production=="other") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz4l-powheg_ntuple.root";
    else assert(0);
  }
  TFile fpowheg(fname);
  TTree *tr_powheg = (TTree*)fpowheg.Get("tr_m4l");
  TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
  tr_powheg->Draw("m4l>>m4l_pow"+binning,"xswgt");
  TH1F *hm = (TH1F*)gDirectory->Get("m4l_pow");
  hm->SetDirectory(0);
  fpowheg.Close();
  hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
  return hm;
}
//----------------------------------------------------------------------------------------
void Draw_em(map<TString,TH1F*> hists, TString proc, TString production, TString mass)
{
  // Draw shapes
  hists["hm_def"]->SetLineWidth(2);
  hists["hm_def"]->Draw();
  hists["hm_lo"]->SetLineColor(kBlue);
  hists["hm_lo"]->SetLineStyle(kDashed);
  hists["hm_lo"]->Draw("same");
  hists["hm_hi"]->SetLineColor(kViolet);
  hists["hm_hi"]->SetLineStyle(kDashed);
  hists["hm_hi"]->Draw("same");
  hists["hpowheg"]->SetLineColor(kRed);
  hists["hpowheg"]->SetLineWidth(2);
  hists["hpowheg"]->Draw("same");
  TLegend leg(0.7,0.7,0.9,0.9);
  leg.AddEntry(hists["hm_def"],"MCFM","L");
  leg.AddEntry(hists["hm_lo"],"MCFM low","L");
  leg.AddEntry(hists["hm_hi"],"MCFM high","L");
  leg.AddEntry(hists["hpowheg"],"POWHEG","L");
  leg.SetFillColor(0);
  leg.Draw();
  // c1->SetLogy();
  c1->SaveAs("./plots/"+proc+"-"+production+"-"+mass+"-shapes.png");

  // Draw ratios
  hists["ratio_def"]->SetLineWidth(2);
  hists["ratio_def"]->Draw();
  hists["ratio_lo"]->SetLineColor(kBlue);
  hists["ratio_lo"]->SetLineStyle(kDashed);
  hists["ratio_lo"]->Draw("same");
  hists["ratio_hi"]->SetLineColor(kRed);
  hists["ratio_hi"]->SetLineStyle(kDashed);
  hists["ratio_hi"]->Draw("same");
  TLegend leg2(0.7,0.7,0.9,0.9);
  leg2.AddEntry(hists["ratio_def"],"MCFM","L");
  leg2.AddEntry(hists["ratio_lo"],"MCFM low","L");
  leg2.AddEntry(hists["ratio_hi"],"MCFM high","L");
  leg2.SetFillColor(0);
  leg2.Draw();
  // c1->SetLogy();
  c1->SaveAs("./plots/"+proc+"-"+production+"-"+mass+"-ratios.png");

}
//----------------------------------------------------------------------------------------
TH1F* makeRatio(TH1F *hm, TH1F *hpowheg)
{
  TH1F *ratio = new TH1F(*hm);
  ratio->Divide(hm,hpowheg);
  for(int ibin=0; ibin<hm->GetNbinsX()+2; ibin++) {
    if(hm->GetBinContent(ibin)==0 || hm->GetBinContent(ibin)==0)
      ratio->SetBinContent(ibin,1);
  }
  return ratio;
}
//----------------------------------------------------------------------------------------
TH1F* makeOppo(TH1F *hm, TH1F *hpowheg)
{
  TH1F *hoppo = new TH1F(*hm);
  for(int ibin=0; ibin<hm->GetNbinsX()+2; ibin++) {
    float powheg = hpowheg->GetBinContent(ibin);
    float mcfm   = hm->GetBinContent(ibin);
    hoppo->SetBinContent(ibin,2*powheg - mcfm);
  }

  return hoppo;
}
//----------------------------------------------------------------------------------------
void writeHists(map<TString,TH1F*> hists, TString process, TString production, TString mass)
{
  TString fname;
  if(process=="higgs")   fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu_weights.root";
  // else if(process=="zz") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz-v11-pu_weights.root";
  else if(process=="zz") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz4l-powheg_weights.root";
  TFile f(fname,"update");
  hists["ratio_def"]->Write("mcfm_def_"+production);
  hists["ratio_lo"]->Write("mcfm_lo_"+production);
  hists["ratio_hi"]->Write("mcfm_hi_"+production);
  f.Close();
}
//----------------------------------------------------------------------------------------
void shapeficator(TString proc, TString production, TString mass)
{
  map<TString,TH1F*> hists;
  hists["hm_def"]  = getMCFM(mass,proc,production,"def");      // mcfm shape
  hists["hm_lo"]   = getMCFM(mass,proc,production,"lo");       // mcfm shape
  hists["hm_hi"]   = getMCFM(mass,proc,production,"hi");       // mcfm shape
  hists["hpowheg"] = getPowheg(mass,proc,production);          // powheg shape
  hists["ratio_def"] = makeRatio(hists["hm_def"],hists["hpowheg"]);
  hists["ratio_lo"]  = makeRatio(hists["hm_lo"], hists["hpowheg"]);
  hists["ratio_hi"]  = makeRatio(hists["hm_hi"], hists["hpowheg"]);
  // cout << mass << "\t" << 3*hists["hpowheg"]->GetRMS() << endl;
  Draw_em(hists,proc,production,mass);
  writeHists(hists,proc,production,mass);
}
//----------------------------------------------------------------------------------------
// main
//----------------------------------------------------------------------------------------
void makeShape()
{
  formatBins();
  shapeficator("zz","gg","115");
  shapeficator("zz","other","115");
  for(unsigned i=0; i<masses.size(); i++) {
    stringstream ss;
    ss << masses[i];
    shapeficator("higgs","gg",ss.str());
  }
}
Revision:	1.3
Committed:	Thu Nov 17 09:45:50 2011 UTC (13 years, 6 months ago) by dkralph
Content type:	text/plain
Branch:	MAIN
CVS Tags:	compiled, synced_FSR_2, synced_FSR, synched2, synched, AN490
Changes since 1.2:	+156 -48 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	dkralph	1.1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
2			#include <TROOT.h> // access to gROOT, entry point to ROOT system
3			#include <TLegend.h>
4			#include <TSystem.h> // interface to OS
5			#include <TFile.h> // file handle class
6			#include <TCanvas.h>
7			#include <TH1F.h>
8			#include <TTree.h> // class to access ntuples
9			#include <TClonesArray.h> // ROOT array class
10			#include <TLorentzVector.h> // 4-vector class
11			#include <TBenchmark.h> // class to track macro running statistics
12			#include <iostream> // standard I/O
13			#include <iomanip> // functions to format standard I/O
14			#include <fstream> // functions for file I/O
15			#include <sstream> // class for parsing strings
16	dkralph	1.3	#include <map>
17	dkralph	1.1
18			#endif
19
20			TCanvas *c1 = new TCanvas("c1","c1");
21	dkralph	1.3	vector<float> masses, widths;
22			map <TString,TString> binstrings;
23			//----------------------------------------------------------------------------------------
24			void formatBins()
25			{
26			binstrings["def"] = "(75,30,500)";
27
28			masses.push_back(115);
29			masses.push_back(120);
30			masses.push_back(130);
31			masses.push_back(140);
32			masses.push_back(150);
33			masses.push_back(160);
34			masses.push_back(170);
35			masses.push_back(180);
36			masses.push_back(190);
37			masses.push_back(200);
38			masses.push_back(230);
39			masses.push_back(250);
40			masses.push_back(300);
41			masses.push_back(350);
42			masses.push_back(400);
43			masses.push_back(500);
44			masses.push_back(600);
45	dkralph	1.1
46	dkralph	1.3	widths.push_back(0.0166566);
47			widths.push_back(0.0185038);
48			widths.push_back(0.0259764);
49			widths.push_back(0.0425892);
50			widths.push_back(0.0920884);
51			widths.push_back(0.438603 );
52			widths.push_back(1.96908 );
53			widths.push_back(3.30333 );
54			widths.push_back(5.40595 );
55			widths.push_back(7.42329 );
56			widths.push_back(14.688 );
57			widths.push_back(20.9409 );
58			widths.push_back(43.5433 );
59			widths.push_back(77.7866 );
60			widths.push_back(141.866 );
61			widths.push_back(266.492 );
62			widths.push_back(379.294 );
63
64			stringstream ss;
65			for(unsigned i=0; i<masses.size(); i++) {
66			float mass = masses[i];
67			float width = widths[i];
68			ss.str("");
69			ss << "(30," << mass-width << "," << mass+width << ")";
70			TString binning(ss.str());
71			ss.str("");
72			ss << mass;
73			binstrings[ss.str()] = binning;
74			}
75			}
76			//----------------------------------------------------------------------------------------
77			TH1F *getMCFM(TString mass, TString process, TString production, TString updown)
78	dkralph	1.1	{
79	dkralph	1.3	TString scale;
80			if(updown=="def") scale = "91_";
81			else if(updown=="lo") scale = "45_";
82			else if(updown=="hi") scale = "182";
83			else assert(0);
84
85			TString fname,wgtStr;
86			if(process=="higgs") {
87			fname = "/temp/dkralph/mcfm/114-"+mass+"-"+updown+"/HZZ_4l_tota_mstw8nl_"+scale+"_"+scale+"_.root";
88			wgtStr = "wt_gg";
89			} else if(process=="zz") {
90			fname = "/temp/dkralph/mcfm/81-115-"+updown+"/ZZlept_tota_mstw8nl_"+scale+"_"+scale+"_.root";
91			if(production=="gg") wgtStr = "wt_gg";
92			else if(production=="other") wgtStr = "wt_gq+wt_qqb";
93			else assert(0);
94			} else assert(0);
95
96			TFile f(fname); assert(f.IsOpen());
97	dkralph	1.1	TTree *tree=0;
98	dkralph	1.3	f.GetObject("h10",tree); assert(tree);
99	dkralph	1.1	TString A,B,C,D,m,varexp;
100			A = "(px3+px4+px5+px6)*(px3+px4+px5+px6)";
101			B = "(py3+py4+py5+py6)*(py3+py4+py5+py6)";
102			C = "(pz3+pz4+pz5+pz6)*(pz3+pz4+pz5+pz6)";
103			D = "(E_3+E_4+E_5+E_6)*(E_3+E_4+E_5+E_6)";
104			m = "sqrt("+D+"-"+A+"-"+B+"-"+C+")";
105	dkralph	1.3	TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
106	dkralph	1.1	varexp= m+">>m4l"+binning;
107	dkralph	1.3	tree->Draw(varexp,wgtStr,"",500000);
108	dkralph	1.1	TH1F hm = (TH1F)gDirectory->Get("m4l");
109			hm->SetDirectory(0);
110			hm->SetTitle("m4l");
111			f.Close();
112			hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
113			return hm;
114			}
115	dkralph	1.2	//----------------------------------------------------------------------------------------
116	dkralph	1.3	TH1F *getPowheg(TString mass, TString process, TString production)
117	dkralph	1.1	{
118			TString fname;
119	dkralph	1.3	if(process=="higgs") {
120			if(production=="gg") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu.root";
121			else assert(0);
122			} else if(process=="zz") {
123			if(production=="gg") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-gg2zz_ntuple.root";
124			else if(production=="other") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz4l-powheg_ntuple.root";
125			else assert(0);
126			}
127	dkralph	1.1	TFile fpowheg(fname);
128			TTree tr_powheg = (TTree)fpowheg.Get("tr_m4l");
129	dkralph	1.3	TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
130	dkralph	1.2	tr_powheg->Draw("m4l>>m4l_pow"+binning,"xswgt");
131	dkralph	1.1	TH1F hm = (TH1F)gDirectory->Get("m4l_pow");
132			hm->SetDirectory(0);
133			fpowheg.Close();
134			hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
135			return hm;
136			}
137	dkralph	1.2	//----------------------------------------------------------------------------------------
138	dkralph	1.3	void Draw_em(map<TString,TH1F*> hists, TString proc, TString production, TString mass)
139	dkralph	1.1	{
140	dkralph	1.3	// Draw shapes
141			hists["hm_def"]->SetLineWidth(2);
142			hists["hm_def"]->Draw();
143			hists["hm_lo"]->SetLineColor(kBlue);
144			hists["hm_lo"]->SetLineStyle(kDashed);
145			hists["hm_lo"]->Draw("same");
146			hists["hm_hi"]->SetLineColor(kViolet);
147			hists["hm_hi"]->SetLineStyle(kDashed);
148			hists["hm_hi"]->Draw("same");
149			hists["hpowheg"]->SetLineColor(kRed);
150			hists["hpowheg"]->SetLineWidth(2);
151			hists["hpowheg"]->Draw("same");
152	dkralph	1.1	TLegend leg(0.7,0.7,0.9,0.9);
153	dkralph	1.3	leg.AddEntry(hists["hm_def"],"MCFM","L");
154			leg.AddEntry(hists["hm_lo"],"MCFM low","L");
155			leg.AddEntry(hists["hm_hi"],"MCFM high","L");
156			leg.AddEntry(hists["hpowheg"],"POWHEG","L");
157	dkralph	1.1	leg.SetFillColor(0);
158	dkralph	1.2	leg.Draw();
159			// c1->SetLogy();
160	dkralph	1.3	c1->SaveAs("./plots/"+proc+"-"+production+"-"+mass+"-shapes.png");
161
162			// Draw ratios
163			hists["ratio_def"]->SetLineWidth(2);
164			hists["ratio_def"]->Draw();
165			hists["ratio_lo"]->SetLineColor(kBlue);
166			hists["ratio_lo"]->SetLineStyle(kDashed);
167			hists["ratio_lo"]->Draw("same");
168			hists["ratio_hi"]->SetLineColor(kRed);
169			hists["ratio_hi"]->SetLineStyle(kDashed);
170			hists["ratio_hi"]->Draw("same");
171			TLegend leg2(0.7,0.7,0.9,0.9);
172			leg2.AddEntry(hists["ratio_def"],"MCFM","L");
173			leg2.AddEntry(hists["ratio_lo"],"MCFM low","L");
174			leg2.AddEntry(hists["ratio_hi"],"MCFM high","L");
175			leg2.SetFillColor(0);
176			leg2.Draw();
177			// c1->SetLogy();
178			c1->SaveAs("./plots/"+proc+"-"+production+"-"+mass+"-ratios.png");
179
180	dkralph	1.1	}
181	dkralph	1.2	//----------------------------------------------------------------------------------------
182	dkralph	1.1	TH1F* makeRatio(TH1F hm, TH1F hpowheg)
183			{
184			TH1F ratio = new TH1F(hm);
185			ratio->Divide(hm,hpowheg);
186	dkralph	1.3	for(int ibin=0; ibin<hm->GetNbinsX()+2; ibin++) {
187			if(hm->GetBinContent(ibin)==0 \|\| hm->GetBinContent(ibin)==0)
188			ratio->SetBinContent(ibin,1);
189			}
190	dkralph	1.1	return ratio;
191			}
192	dkralph	1.2	//----------------------------------------------------------------------------------------
193	dkralph	1.1	TH1F* makeOppo(TH1F hm, TH1F hpowheg)
194			{
195			TH1F hoppo = new TH1F(hm);
196			for(int ibin=0; ibin<hm->GetNbinsX()+2; ibin++) {
197			float powheg = hpowheg->GetBinContent(ibin);
198			float mcfm = hm->GetBinContent(ibin);
199			hoppo->SetBinContent(ibin,2*powheg - mcfm);
200			}
201
202			return hoppo;
203			}
204	dkralph	1.2	//----------------------------------------------------------------------------------------
205	dkralph	1.3	void writeHists(map<TString,TH1F*> hists, TString process, TString production, TString mass)
206	dkralph	1.1	{
207			TString fname;
208			if(process=="higgs") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu_weights.root";
209	dkralph	1.2	// else if(process=="zz") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz-v11-pu_weights.root";
210			else if(process=="zz") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz4l-powheg_weights.root";
211	dkralph	1.1	TFile f(fname,"update");
212	dkralph	1.3	hists["ratio_def"]->Write("mcfm_def_"+production);
213			hists["ratio_lo"]->Write("mcfm_lo_"+production);
214			hists["ratio_hi"]->Write("mcfm_hi_"+production);
215	dkralph	1.1	f.Close();
216			}
217			//----------------------------------------------------------------------------------------
218	dkralph	1.3	void shapeficator(TString proc, TString production, TString mass)
219			{
220			map<TString,TH1F*> hists;
221			hists["hm_def"] = getMCFM(mass,proc,production,"def"); // mcfm shape
222			hists["hm_lo"] = getMCFM(mass,proc,production,"lo"); // mcfm shape
223			hists["hm_hi"] = getMCFM(mass,proc,production,"hi"); // mcfm shape
224			hists["hpowheg"] = getPowheg(mass,proc,production); // powheg shape
225			hists["ratio_def"] = makeRatio(hists["hm_def"],hists["hpowheg"]);
226			hists["ratio_lo"] = makeRatio(hists["hm_lo"], hists["hpowheg"]);
227			hists["ratio_hi"] = makeRatio(hists["hm_hi"], hists["hpowheg"]);
228			// cout << mass << "\t" << 3*hists["hpowheg"]->GetRMS() << endl;
229			Draw_em(hists,proc,production,mass);
230			writeHists(hists,proc,production,mass);
231			}
232			//----------------------------------------------------------------------------------------
233	dkralph	1.1	// main
234			//----------------------------------------------------------------------------------------
235			void makeShape()
236			{
237	dkralph	1.3	formatBins();
238			shapeficator("zz","gg","115");
239			shapeficator("zz","other","115");
240			for(unsigned i=0; i<masses.size(); i++) {
241			stringstream ss;
242			ss << masses[i];
243			shapeficator("higgs","gg",ss.str());
244			}
245	dkralph	1.1	}