ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitHzz4l/Shapes/makeShape.C

Comparing UserCode/MitHzz4l/Shapes/makeShape.C (file contents):
Revision 1.2 by dkralph, Tue Nov 15 07:46:59 2011 UTC vs.
Revision 1.3 by dkralph, Thu Nov 17 09:45:50 2011 UTC

#	Line 13 | Line 13
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	+	#include <map>
17
18		#endif
19
20		TCanvas *c1 = new TCanvas("c1","c1");
21	<	TString binning("(100,0,500)");
21	>	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
46	<	TH1F *getMCFM(TString mass, TString process)
46	>	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		{
79	<	TString fname;
80	<	if(process=="higgs")
81	<	fname = "/temp/dkralph/mcfm/114-"+mass+"/HZZ_4l_virt_mstw8nl_91__91__virtggonly.root";
82	<	else if(process=="zz")
83	<	fname = "/temp/dkralph/mcfm/81-115/ZZlept_tota_mstw8nl_91__91__virtggonly.root";
84	<	// fname = "/temp/dkralph/ZZlept_virt_mstw8nl_91__91__virtggonly.root";
85	<	TFile f(fname);
79	>	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		TTree *tree=0;
98	<	f.GetObject("h10",tree);
98	>	f.GetObject("h10",tree); assert(tree);
99		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	+	TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
106		varexp= m+">>m4l"+binning;
107	<	tree->Draw(varexp,"wt_ALL");
107	>	tree->Draw(varexp,wgtStr,"",500000);
108		TH1F *hm = (TH1F*)gDirectory->Get("m4l");
109		hm->SetDirectory(0);
110		hm->SetTitle("m4l");
#	Line 46 | Line 113 | TH1F *getMCFM(TString mass, TString proc
113		return hm;
114		}
115		//----------------------------------------------------------------------------------------
116	<	TH1F *getPowheg(TString mass, TString process)
116	>	TH1F *getPowheg(TString mass, TString process, TString production)
117		{
118		TString fname;
119	<	if(process=="higgs")
120	<	fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu.root";
121	<	if(process=="zz")
122	<	// fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz-v11-pu_ntuple.root";
123	<	fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz4l-powheg_ntuple.root";
124	<
119	>	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		TFile fpowheg(fname);
128		TTree *tr_powheg = (TTree*)fpowheg.Get("tr_m4l");
129	+	TString binning = process=="zz" ? binstrings["def"] : binstrings[mass];
130		tr_powheg->Draw("m4l>>m4l_pow"+binning,"xswgt");
131		TH1F *hm = (TH1F*)gDirectory->Get("m4l_pow");
132		hm->SetDirectory(0);
#	Line 65 | Line 135 | TH1F *getPowheg(TString mass, TString pr
135		return hm;
136		}
137		//----------------------------------------------------------------------------------------
138	<	void Draw_em(TH1F *hm, TH1F *hpowheg, TH1F *hratio, TH1F *hoppo, TH1F *hopporatio)
138	>	void Draw_em(map<TString,TH1F*> hists, TString proc, TString production, TString mass)
139		{
140	<	hm->Draw();
141	<	hpowheg->SetLineColor(kRed);
142	<	hpowheg->Draw("same");
143	<	hratio->SetLineColor(kBlue);
144	<	// hratio->Scale(1./hratio->Integral(0,hratio->GetNbinsX()));  // don't scale it!
145	<	// hratio->Draw("same");
146	<	hoppo->SetLineColor(kBlue);
147	<	hoppo->Scale(1./hoppo->Integral(0,hoppo->GetNbinsX()));
148	<	// hoppo->Draw("same");
149	<	hopporatio->SetLineColor(kViolet);
150	<	// hopporatio->Scale(1./hopporatio->Integral(0,hopporatio->GetNbinsX())); // don't scale it!
151	<	// hopporatio->Draw("same");
140	>	// 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		TLegend leg(0.7,0.7,0.9,0.9);
153	<	leg.AddEntry(hm,"MCFM","L");
154	<	leg.AddEntry(hpowheg,"POWHEG","L");
155	<	leg.AddEntry(hratio,"ratio","L");
156	<	leg.AddEntry(hoppo,"opp. mcfm","L");
87	<	leg.AddEntry(hopporatio,"opp. ratio","L");
153	>	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		leg.SetFillColor(0);
158		leg.Draw();
159		// c1->SetLogy();
160	<	c1->SaveAs("~/public_html/fo3o.png");
160	>	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		}
181		//----------------------------------------------------------------------------------------
182		TH1F* makeRatio(TH1F *hm, TH1F *hpowheg)
183		{
184		TH1F *ratio = new TH1F(*hm);
185		ratio->Divide(hm,hpowheg);
186	+	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		return ratio;
191		}
192		//----------------------------------------------------------------------------------------
#	Line 110 | Line 202 | TH1F* makeOppo(TH1F *hm, TH1F *hpowheg)
202		return hoppo;
203		}
204		//----------------------------------------------------------------------------------------
205	<	void writeHists(TH1F *hratio, TH1F *hopporatio, TString process, TString mass)
205	>	void writeHists(map<TString,TH1F*> hists, TString process, TString production, TString mass)
206		{
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		// 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		TFile f(fname,"update");
212	<	hratio->Write("low");
213	<	hopporatio->Write("high");
212	>	hists["ratio_def"]->Write("mcfm_def_"+production);
213	>	hists["ratio_lo"]->Write("mcfm_lo_"+production);
214	>	hists["ratio_hi"]->Write("mcfm_hi_"+production);
215		f.Close();
216		}
217		//----------------------------------------------------------------------------------------
218	+	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		// main
234		//----------------------------------------------------------------------------------------
235		void makeShape()
236		{
237	<	TString mass("115");
238	<	TH1F *hm = getMCFM(mass,"zz");                 // mcfm shape
239	<	TH1F *hpowheg = getPowheg(mass,"zz");          // powheg shape
240	<	TH1F *hratio = makeRatio(hm,hpowheg);          // mcfm/powheg
241	<	TH1F *hoppo  = makeOppo(hm,hpowheg);           // "opposite" shape to mcfm
242	<	TH1F *hopporatio = makeRatio(hoppo,hpowheg);   // "opposite to mcfm"/powheg
243	<	Draw_em(hm,hpowheg,hratio,hoppo,hopporatio);
244	<	writeHists(hratio,hopporatio,"zz",mass);
237	>	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		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines