ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Shapes/macros/makeShape.C
Revision: 1.2
Committed: Fri Dec 2 17:09:03 2011 UTC (13 years, 5 months ago) by dkralph
Content type: text/plain
Branch: MAIN
CVS Tags: compiled, synced_FSR_2, synced_FSR, synched2, synched, AN490, HEAD
Changes since 1.1: +32 -26 lines
Log Message:
*** empty log message ***

File Contents

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