ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Shapes/makeShape.C
Revision: 1.1
Committed: Mon Oct 31 13:51:42 2011 UTC (13 years, 6 months ago) by dkralph
Content type: text/plain
Branch: MAIN
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 dkralph 1.1 // #define debug
2     //================================================================================================
3     //
4     // Measure fake rate using muon triggered sample
5     //
6     //________________________________________________________________________________________________
7    
8     #if !defined(__CINT__) || defined(__MAKECINT__)
9     #include <TROOT.h> // access to gROOT, entry point to ROOT system
10     #include <TLegend.h>
11     #include <TSystem.h> // interface to OS
12     #include <TFile.h> // file handle class
13     #include <TCanvas.h>
14     #include <TH1F.h>
15     #include <TTree.h> // class to access ntuples
16     #include <TClonesArray.h> // ROOT array class
17     #include <TLorentzVector.h> // 4-vector class
18     #include <TBenchmark.h> // class to track macro running statistics
19     #include <iostream> // standard I/O
20     #include <iomanip> // functions to format standard I/O
21     #include <fstream> // functions for file I/O
22     #include <sstream> // class for parsing strings
23    
24     #endif
25    
26     TCanvas *c1 = new TCanvas("c1","c1");
27     TString binning("(100,90,500)");
28    
29     TH1F *getMCFM(TString mass, TString process)
30     {
31     TString fname;
32     if(process=="higgs")
33     fname = "/temp/dkralph/mcfm/114-"+mass+"/HZZ_4l_virt_mstw8nl_91__91__virtggonly.root";
34     else if(process=="zz")
35     fname = "/temp/dkralph/mcfm/81-115/ZZlept_tota_mstw8nl_91__91__virtggonly.root";
36     // fname = "/temp/dkralph/ZZlept_virt_mstw8nl_91__91__virtggonly.root";
37     TFile f(fname);
38     TTree *tree=0;
39     f.GetObject("h10",tree);
40     TString A,B,C,D,m,varexp;
41     A = "(px3+px4+px5+px6)*(px3+px4+px5+px6)";
42     B = "(py3+py4+py5+py6)*(py3+py4+py5+py6)";
43     C = "(pz3+pz4+pz5+pz6)*(pz3+pz4+pz5+pz6)";
44     D = "(E_3+E_4+E_5+E_6)*(E_3+E_4+E_5+E_6)";
45     m = "sqrt("+D+"-"+A+"-"+B+"-"+C+")";
46     varexp= m+">>m4l"+binning;
47     tree->Draw(varexp,"wt_ALL");
48     TH1F *hm = (TH1F*)gDirectory->Get("m4l");
49     hm->SetDirectory(0);
50     hm->SetTitle("m4l");
51     f.Close();
52     hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
53     return hm;
54     }
55    
56     TH1F *getPowheg(TString mass, TString process)
57     {
58     TString fname;
59     if(process=="higgs")
60     fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu.root";
61     if(process=="zz")
62     fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz-v11-pu_ntuple.root";
63    
64     TFile fpowheg(fname);
65     TTree *tr_powheg = (TTree*)fpowheg.Get("tr_m4l");
66     tr_powheg->Draw("m4l>>m4l_pow"+binning);
67     TH1F *hm = (TH1F*)gDirectory->Get("m4l_pow");
68     hm->SetDirectory(0);
69     fpowheg.Close();
70     hm->Scale(1./hm->Integral(0,hm->GetNbinsX()));
71     return hm;
72     }
73     void Draw_em(TH1F *hm, TH1F *hpowheg, TH1F *hratio, TH1F *hoppo, TH1F *hopporatio)
74     {
75     hm->Draw();
76     hpowheg->SetLineColor(kRed);
77     hpowheg->Draw("same");
78     hratio->SetLineColor(kBlue);
79     // hratio->Scale(1./hratio->Integral(0,hratio->GetNbinsX())); // don't scale it!
80     hratio->Draw("same");
81     hoppo->SetLineColor(kBlack);
82     hoppo->Scale(1./hoppo->Integral(0,hoppo->GetNbinsX()));
83     hoppo->Draw("same");
84     hopporatio->SetLineColor(kViolet);
85     // hopporatio->Scale(1./hopporatio->Integral(0,hopporatio->GetNbinsX())); // don't scale it!
86     hopporatio->Draw("same");
87     TLegend leg(0.7,0.7,0.9,0.9);
88     leg.AddEntry(hm,"MCFM","L");
89     leg.AddEntry(hpowheg,"POWHEG","L");
90     leg.AddEntry(hratio,"ratio","L");
91     leg.AddEntry(hoppo,"opp. mcfm","L");
92     leg.AddEntry(hopporatio,"opp. ratio","L");
93     leg.SetFillColor(0);
94     // leg.Draw();
95     c1->SetLogy();
96     c1->SaveAs("~/public_html/fo3o.png");
97     }
98     TH1F* makeRatio(TH1F *hm, TH1F *hpowheg)
99     {
100     TH1F *ratio = new TH1F(*hm);
101     ratio->Divide(hm,hpowheg);
102     return ratio;
103     }
104     TH1F* makeOppo(TH1F *hm, TH1F *hpowheg)
105     {
106     TH1F *hoppo = new TH1F(*hm);
107     for(int ibin=0; ibin<hm->GetNbinsX()+2; ibin++) {
108     float powheg = hpowheg->GetBinContent(ibin);
109     float mcfm = hm->GetBinContent(ibin);
110     hoppo->SetBinContent(ibin,2*powheg - mcfm);
111     }
112    
113     return hoppo;
114     }
115     void writeHists(TH1F *hratio, TH1F *hopporatio, TString process, TString mass)
116     {
117     TString fname;
118     if(process=="higgs") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-h"+mass+"zz4l-gf-v1g1-pu_weights.root";
119     else if(process=="zz") fname = "~/cms/cmssw/025/CMSSW_4_4_1/src/MitHzz4l/data/m4lshapes/s11-zz-v11-pu_weights.root";
120     TFile f(fname,"update");
121     hratio->Write("low");
122     hopporatio->Write("high");
123     f.Close();
124     }
125     //----------------------------------------------------------------------------------------
126     // main
127     //----------------------------------------------------------------------------------------
128     void makeShape()
129     {
130     TString mass("115");
131     TH1F *hm = getMCFM(mass,"zz"); // mcfm shape
132     TH1F *hpowheg = getPowheg(mass,"zz"); // powheg shape
133     TH1F *hratio = makeRatio(hm,hpowheg); // mcfm/powheg
134     TH1F *hoppo = makeOppo(hm,hpowheg); // "opposite" shape to mcfm
135     TH1F *hopporatio = makeRatio(hoppo,hpowheg); // "opposite to mcfm"/powheg
136     Draw_em(hm,hpowheg,hratio,hoppo,hopporatio);
137     writeHists(hratio,hopporatio,"zz",mass);
138     }