ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ResultModule.C
Revision: 1.3
Committed: Tue Jul 12 09:00:15 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.2: +19 -7 lines
Log Message:
Updated result module, added poisson statistics as well as  stat & syst errors

File Contents

# User Rev Content
1 buchmann 1.1 #include <iostream>
2     #include <vector>
3     #include <sys/stat.h>
4    
5     #include <TCut.h>
6     #include <TROOT.h>
7     #include <TCanvas.h>
8     #include <TMath.h>
9     #include <TColor.h>
10     #include <TPaveText.h>
11     #include <TRandom.h>
12     #include <TH1.h>
13     #include <TH2.h>
14     #include <TF1.h>
15     #include <TSQLResult.h>
16    
17     //#include "TTbar_stuff.C"
18     using namespace std;
19    
20     using namespace PlottingSetup;
21    
22 buchmann 1.2 string newjzbexpression(string oldexpression,float shift) {
23     stringstream ss;
24     if(shift>0) ss<<oldexpression<<"+"<<shift;
25     if(shift<0) ss<<oldexpression<<shift;
26     if(shift==0) ss<<oldexpression;
27     return ss.str();
28     }
29 buchmann 1.1
30 buchmann 1.2 void fill_result_histos(float &zossfp, float &zosofp, float &zossfn, float &zosofn, float &sbossfp, float &sbosofp, float &sbossfn, float &sbosofn,string datajzb,float cut, int mcordata, float &result) {
31 buchmann 1.1 string xlabel="JZB (GeV) -- for algoritm internal use only!";
32     TH1F *ZOSSFP = allsamples.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
33     TH1F *ZOSOFP = allsamples.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
34     TH1F *ZOSSFN = allsamples.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccut,mcordata,luminosity);
35     TH1F *ZOSOFN = allsamples.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccut,mcordata,luminosity);
36    
37     TH1F *SBOSSFP = allsamples.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
38     TH1F *SBOSOFP = allsamples.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
39     TH1F *SBOSSFN = allsamples.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
40     TH1F *SBOSOFN = allsamples.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccut&&sidebandcut,mcordata,luminosity);
41 buchmann 1.2
42     zossfp=ZOSSFP->Integral();
43     zosofp=ZOSOFP->Integral();
44     zossfn=ZOSSFN->Integral();
45     zosofn=ZOSOFN->Integral();
46     sbossfp=SBOSSFP->Integral();
47     sbosofp=SBOSOFP->Integral();
48     sbossfn=SBOSSFN->Integral();
49     sbosofn=SBOSOFN->Integral();
50    
51     result = zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn);
52     delete ZOSSFP;
53     delete ZOSOFP;
54     delete ZOSSFN;
55     delete ZOSOFN;
56    
57     delete SBOSSFP;
58     delete SBOSOFP;
59     delete SBOSSFN;
60     delete SBOSOFN;
61    
62     }
63    
64    
65 buchmann 1.3 void get_result_above_one_fixed_jzb_value(float cut ,string mcjzb,string datajzb, int mcordata,float jzbpeakerrorMC, float jzbpeakerrorData, TCanvas *rescan, bool chatty=false, bool dopoisson=false) {
66 buchmann 1.2 rescan->cd();
67     if(mcordata!=mc) cout << "Crunching the numbers for JZB>" << cut << endl;
68    
69     float zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,result;
70     if(mcordata==mc) fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,mcordata,result);
71     if(mcordata==data) fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,mcordata,result);
72    
73    
74     float ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn, ppresult;
75     if(mcordata==mc) fill_result_histos(ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,mcordata,ppresult);
76     if(mcordata==data) fill_result_histos(ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(datajzb,jzbpeakerrorData),cut,mcordata,ppresult);
77    
78     float pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn, pnresult;
79     if(mcordata==mc) fill_result_histos(pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,mcordata,pnresult);
80     if(mcordata==data) fill_result_histos(pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(datajzb,-jzbpeakerrorData),cut,mcordata,pnresult);
81 buchmann 1.1
82     float syserr=0;
83     float peakerr=0;
84     float staterr=0;
85 buchmann 1.3 float poissonstaterrup=-999;
86     float poissonstaterrdown=-999;
87     if(fabs(result-pnresult)>fabs(result-ppresult)) peakerr=fabs(result-pnresult); else peakerr=fabs(result-ppresult);
88 buchmann 1.1
89 buchmann 1.3 float zjetsestimateuncert=0.5;
90     float emuncert=0.5;
91     float emsidebanduncert=0.5;
92     float eemmsidebanduncert=0.5;
93     syserr = (zjetsestimateuncert*zossfn)*(zjetsestimateuncert*zossfn);//first term
94     syserr+= ((zosofp)*(zosofp) + (zosofn)*(zosofn))*(1.0/9)*emuncert*emuncert;//sys err from emu method
95     syserr+= ((sbossfp)*(sbossfp)+(sbossfn)*(sbossfn))*(1.0/9)*eemmsidebanduncert*eemmsidebanduncert; // sys err from eemm sidebands
96     syserr+= ((sbosofp)*(sbosofp)+(sbosofn)*(sbosofn))*(1.0/9)*emsidebanduncert*emsidebanduncert; // sys err from emu sidebands
97     syserr=TMath::Sqrt(syserr);
98 buchmann 1.2
99 buchmann 1.3 staterr=TMath::Sqrt(zossfn + (1.0/9)*(zosofp+zosofn)+ (1.0/9)*(sbossfp+sbossfn)+ (1.0/9)*(sbosofp+sbosofn));
100     if(dopoisson) advanced_poisson(zossfn,zosofp,zosofn,sbossfp,sbossfn,sbosofp,sbosofn,poissonstaterrdown,poissonstaterrup);
101 buchmann 1.2
102 buchmann 1.1 if(mcordata==mc) cout << " MC :: ";
103     else cout << " ";
104 buchmann 1.2 cout << "Observed : " << zossfp << endl;
105 buchmann 1.1 if(mcordata==mc) cout << " MC :: ";
106     else cout << " ";
107 buchmann 1.2 cout << "Predicted: " << zossfn << " + (1/3)*(" << zosofp << "-" << zosofn<<") [e&mu]+ (1/3)*(" << sbossfp << "-" << sbossfn<<") [SF,SB]+ (1/3)*(" << sbosofp << "-" << sbosofn<<") [OF,SB] = ";
108     cout << zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn)<< " +/- " << peakerr << " (peak ) " << " +/- " << syserr << " (sys) +/- " << staterr << " (stat)" << endl;
109 buchmann 1.3 if(dopoisson) cout << " stat error with Poisson : +" << poissonstaterrup << " - " << poissonstaterrdown << endl;
110 buchmann 1.1
111     if(chatty) {
112 buchmann 1.2 cout << " Pred(ZJets ) \t " << zossfn << endl;
113     cout << " Pred(e&mu;]) \t " << zosofp << "-" << zosofn << " = " << zosofp-zosofn<<endl;
114     cout << " Pred(ossf,sb]) \t " << sbossfp << "-" << sbossfn<<" = "<<sbossfp-sbossfn<<endl;
115     cout << " Pred(osof,sb]) \t " << sbosofp << "-" << sbosofn<<" = "<<sbosofp-sbosofn<<endl;
116 buchmann 1.1 }
117     }
118    
119    
120 buchmann 1.3 void get_result(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbcuts, bool chatty=false, bool dopoisson=false) {
121 buchmann 1.1 TCanvas *rescan = new TCanvas("rescan","Result Canvas");
122     //void look_at_sidebands(string mcjzb, string datajzb) {
123     for(int icut=0;icut<jzbcuts.size();icut++) {
124 buchmann 1.3 get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,data,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
125     get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,mc,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
126 buchmann 1.1 }
127     delete rescan;
128     }