ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ResultModule.C
Revision: 1.14
Committed: Wed Sep 14 06:59:14 2011 UTC (13 years, 7 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: Honeypot, cbaf_2p1ifb
Changes since 1.13: +2 -0 lines
Log Message:
Offpeak warnings added

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.11 void fill_result_histos(float &zossfp, float &zossfperr, float &zosofp, float &zossfn, float &zosofn, float &sbossfp, float &sbosofp, float &sbossfn, float &sbosofn,string datajzb,float cut, float cuthigh, int mcordata, float &result, vector<int> sel, samplecollection &sampleC, string addcut="") {
23 buchmann 1.14 if(!RestrictToMassPeak) write_warning(__FUNCTION__,"Watch out, once we go offpeak this function needs rewriting!");
24 buchmann 1.5 string xlabel="JZB [GeV] -- for algoritm internal use only!";
25     bool dosignal=false;
26     if(mcordata==mcwithsignal) {
27     dosignal=true;
28     mcordata=0;
29     }
30 buchmann 1.8 TCut basiccutplus;
31     if(addcut=="") basiccutplus=basiccut;
32     else basiccutplus=basiccut&&addcut.c_str();
33 buchmann 1.6 TH1F *ZOSSFP;
34     TH1F *ZOSOFP;
35     TH1F *ZOSSFN;
36     TH1F *ZOSOFN;
37    
38     TH1F *SBOSSFP;
39     TH1F *SBOSOFP;
40     TH1F *SBOSSFN;
41     TH1F *SBOSOFN;
42    
43     if(mcordata==mc||mcordata==data||mcordata==mcwithsignal) {
44 buchmann 1.11 ZOSSFP = sampleC.Draw("ZOSSFP",datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
45     ZOSOFP = sampleC.Draw("ZOSOFP",datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
46     ZOSSFN = sampleC.Draw("ZOSSFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
47     ZOSOFN = sampleC.Draw("ZOSOFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
48    
49     SBOSSFP = sampleC.Draw("SBOSSFP",datajzb,1,cut,cuthigh, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
50     SBOSOFP = sampleC.Draw("SBOSOFP",datajzb,1,cut,cuthigh, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
51     SBOSSFN = sampleC.Draw("SBOSSFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
52     SBOSOFN = sampleC.Draw("SBOSOFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
53 buchmann 1.6 } else {
54     //doing signal only!
55 buchmann 1.11 ZOSSFP = sampleC.Draw("ZOSSFP",datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
56     ZOSOFP = sampleC.Draw("ZOSOFP",datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
57     ZOSSFN = sampleC.Draw("ZOSSFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
58     ZOSOFN = sampleC.Draw("ZOSOFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
59    
60     SBOSSFP = sampleC.Draw("SBOSSFP",datajzb,1,cut,cuthigh, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
61     SBOSOFP = sampleC.Draw("SBOSOFP",datajzb,1,cut,cuthigh, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
62     SBOSSFN = sampleC.Draw("SBOSSFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
63     SBOSOFN = sampleC.Draw("SBOSOFN","-"+datajzb,1,cut,cuthigh, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
64 buchmann 1.6 }
65    
66 fronga 1.9 double err; // UGH!
67 buchmann 1.10 zossfp=IntegralAndError(ZOSSFP,1,ZOSSFP->GetNbinsX(),err,"");//making this compatible with my computer which has an outdated version of root.
68 fronga 1.9 zossfperr = err;
69 buchmann 1.5 zosofp=ZOSOFP->Integral();
70     zossfn=ZOSSFN->Integral();
71     zosofn=ZOSOFN->Integral();
72    
73     sbossfp=SBOSSFP->Integral();
74     sbosofp=SBOSOFP->Integral();
75     sbossfn=SBOSSFN->Integral();
76     sbosofn=SBOSOFN->Integral();
77 buchmann 1.2
78     delete ZOSSFP;
79     delete ZOSOFP;
80     delete ZOSSFN;
81     delete ZOSOFN;
82    
83     delete SBOSSFP;
84     delete SBOSOFP;
85     delete SBOSSFN;
86     delete SBOSOFN;
87 buchmann 1.6
88 fronga 1.9 result = zossfn + (1.0/3.)*(zosofp-zosofn)+ (1.0/3.)*(sbossfp-sbossfn)+ (1.0/3.)*(sbosofp-sbosofn);
89 buchmann 1.6 }
90    
91 buchmann 1.11 void fill_result_histos(float &zossfp, float &zossfperr, float &zosofp, float &zossfn, float &zosofn, float &sbossfp, float &sbosofp, float &sbossfn, float &sbosofn,string datajzb,float cut, float cuthigh, int mcordata, float &result, string addcut="") {
92 buchmann 1.6 vector<int> emptyvector;
93 buchmann 1.11 fill_result_histos(zossfp, zossfperr,zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,cuthigh,mcordata,result,emptyvector,allsamples,addcut);
94 buchmann 1.2 }
95    
96 buchmann 1.11 /*
97 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) {
98 buchmann 1.2 rescan->cd();
99 fronga 1.9 if(mcordata==data) dout << "***\nCrunching numbers for JZB>" << cut << " : " << endl;
100 buchmann 1.2
101 fronga 1.9 float zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,result;
102     if(mcordata==mc) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,mcordata,result);
103     if(mcordata==data) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,mcordata,result);
104     if(mcordata==mcwithsignal) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,mcordata,result);
105 buchmann 1.2
106 fronga 1.9 float zossfpee, zossfpeeerr, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,resultee;
107     float zossfpmm, zossfpmmerr, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,resultmm;
108 buchmann 1.8 if(mcordata==data) {
109 fronga 1.9 fill_result_histos(zossfpee, zossfpeeerr, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,datajzb,cut,mcordata,resultee,"id1==0");
110     fill_result_histos(zossfpmm, zossfpmmerr, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,datajzb,cut,mcordata,resultmm,"id1==1");
111 buchmann 1.8 }
112    
113 fronga 1.9 float ppzossfp, ppzossfperr, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn, ppresult;
114     if(mcordata==mc) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,mcordata,ppresult);
115     if(mcordata==data) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(datajzb,jzbpeakerrorData),cut,mcordata,ppresult);
116     if(mcordata==mcwithsignal) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,mcordata,ppresult);
117    
118     float pnzossfp, pnzossfperr, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn, pnresult;
119     if(mcordata==mc) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,mcordata,pnresult);
120     if(mcordata==data) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(datajzb,-jzbpeakerrorData),cut,mcordata,pnresult);
121     if(mcordata==mcwithsignal) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,mcordata,pnresult);
122 buchmann 1.1
123     float syserr=0;
124     float peakerr=0;
125     float staterr=0;
126 buchmann 1.3 float poissonstaterrup=-999;
127     float poissonstaterrdown=-999;
128     if(fabs(result-pnresult)>fabs(result-ppresult)) peakerr=fabs(result-pnresult); else peakerr=fabs(result-ppresult);
129 buchmann 1.1
130 buchmann 1.3 float zjetsestimateuncert=0.5;
131     float emuncert=0.5;
132     float emsidebanduncert=0.5;
133     float eemmsidebanduncert=0.5;
134     syserr = (zjetsestimateuncert*zossfn)*(zjetsestimateuncert*zossfn);//first term
135     syserr+= ((zosofp)*(zosofp) + (zosofn)*(zosofn))*(1.0/9)*emuncert*emuncert;//sys err from emu method
136     syserr+= ((sbossfp)*(sbossfp)+(sbossfn)*(sbossfn))*(1.0/9)*eemmsidebanduncert*eemmsidebanduncert; // sys err from eemm sidebands
137     syserr+= ((sbosofp)*(sbosofp)+(sbosofn)*(sbosofn))*(1.0/9)*emsidebanduncert*emsidebanduncert; // sys err from emu sidebands
138     syserr=TMath::Sqrt(syserr);
139 buchmann 1.2
140 buchmann 1.3 staterr=TMath::Sqrt(zossfn + (1.0/9)*(zosofp+zosofn)+ (1.0/9)*(sbossfp+sbossfn)+ (1.0/9)*(sbosofp+sbosofn));
141     if(dopoisson) advanced_poisson(zossfn,zosofp,zosofn,sbossfp,sbossfn,sbosofp,sbosofn,poissonstaterrdown,poissonstaterrup);
142 buchmann 1.2
143 buchmann 1.8 float e_to_emu=0.5;
144     float m_to_emu=0.5;
145    
146 buchmann 1.7 if(mcordata==mc) dout << " MC :: ";
147     if(mcordata==mcwithsignal) dout << " MC with S :: ";
148     if(mcordata==data) dout << " ";
149 fronga 1.9 dout << "Observed : " << zossfp << "+-" << zossfperr << endl;
150 buchmann 1.8 if(mcordata==data) dout << " Composition: " << zossfpee << " (ee), " << zossfpmm << " (mm) " << endl;
151    
152 buchmann 1.7 if(mcordata==mc) dout << " MC :: ";
153     if(mcordata==mcwithsignal) dout << " MC with S :: ";
154     if(mcordata==data) dout << " ";
155     dout << "Predicted: " << zossfn << " + (1/3)*(" << zosofp << "-" << zosofn<<") [e&mu]+ (1/3)*(" << sbossfp << "-" << sbossfn<<") [SF,SB]+ (1/3)*(" << sbosofp << "-" << sbosofn<<") [OF,SB] = ";
156     dout << zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn)<< " +/- " << peakerr << " (peak ) " << " +/- " << syserr << " (sys) +/- " << staterr << " (stat)" << endl;
157     if(dopoisson) dout << " stat error with Poisson : +" << poissonstaterrup << " - " << poissonstaterrdown << endl;
158 buchmann 1.1
159 buchmann 1.8 if(mcordata==data) {
160     dout << " Composition: (ee)" << zossfnee << " + (1/3)*(" << e_to_emu*zosofp << "-" << e_to_emu*zosofn<<") ["<<e_to_emu<<"*e&mu]+ (1/3)*(" << sbossfpee << "-" << sbossfnee<<") [SF,SB]+ (1/3)*(" << sbosofp*e_to_emu << "-" << sbosofn*e_to_emu<<") ["<<e_to_emu<<"*OF,SB] = ";
161     dout << zossfnee + (1.0/3)*(e_to_emu*zosofp-e_to_emu*zosofn)+ (1.0/3)*(sbossfpee-sbossfnee) + (1.0/3)*(sbosofp*e_to_emu - sbosofn*e_to_emu) << endl;
162     dout << " (mm)" << zossfnmm << " + (1/3)*(" << m_to_emu*zosofp << "-" << m_to_emu*zosofn<<") ["<<m_to_emu<<"*e&mu]+ (1/3)*(" << sbossfpmm << "-" << sbossfnmm<<") [SF,SB]+ (1/3)*(" << sbosofp*m_to_emu << "-" << sbosofn*m_to_emu<<") ["<<m_to_emu<<"*OF,SB] = ";
163     dout << zossfnmm + (1.0/3)*(m_to_emu*zosofp - m_to_emu*zosofn) + (1.0/3)*(sbossfpmm -sbossfnmm)+ (1.0/3)*(sbosofp*m_to_emu-sbosofn*e_to_emu) << endl;
164     }
165    
166 buchmann 1.1 if(chatty) {
167 buchmann 1.7 dout << " Pred(ZJets ) \t " << zossfn << endl;
168     dout << " Pred(e&mu;]) \t " << zosofp << "-" << zosofn << " = " << zosofp-zosofn<<endl;
169     dout << " Pred(ossf,sb]) \t " << sbossfp << "-" << sbossfn<<" = "<<sbossfp-sbossfn<<endl;
170     dout << " Pred(osof,sb]) \t " << sbosofp << "-" << sbosofn<<" = "<<sbosofp-sbosofn<<endl;
171 buchmann 1.1 }
172 buchmann 1.6
173     if(mcordata==data) {
174     //store the result!
175     Nobs.push_back(zossfp);
176     Npred.push_back(zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn));
177     float totprederr=0;
178     totprederr+=peakerr*peakerr;
179     totprederr+=syserr*syserr;
180     totprederr+=staterr*staterr;
181     totprederr=TMath::Sqrt(totprederr);
182     Nprederr.push_back(totprederr);
183     }
184 buchmann 1.1 }
185 buchmann 1.11 */
186     vector<float> get_result_between_two_fixed_jzb_values(float cut , float cuthigh, string mcjzb,string datajzb, int mcordata,float jzbpeakerrorMC, float jzbpeakerrorData, TCanvas *rescan, bool chatty=false, bool dopoisson=false, bool writeanything=true) {
187 buchmann 1.14 if(!RestrictToMassPeak) write_warning(__FUNCTION__,"Watch out, once we go offpeak this function needs rewriting!");
188 buchmann 1.11 /*return vector of floats
189     [0] Bpred [1] Bpred uncert [2] Observed [3] Observed error
190     // if we use this for the ratio plot we don't want to see any of the results (hence writeanything=false)
191     */
192     rescan->cd();
193     if(writeanything&&mcordata==data) dout << "***\nCrunching numbers for JZB>" << cut << " : " << endl;
194    
195     float zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,result;
196     if(mcordata==mc) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,cuthigh,mcordata,result);
197     if(mcordata==data) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,cuthigh,mcordata,result);
198     if(mcordata==mcwithsignal) fill_result_histos(zossfp, zossfperr, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,cuthigh,mcordata,result);
199    
200     float zossfpee, zossfpeeerr, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,resultee;
201     float zossfpmm, zossfpmmerr, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,resultmm;
202     if(mcordata==data) {
203     fill_result_histos(zossfpee, zossfpeeerr, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,datajzb,cut,cuthigh,mcordata,resultee,"id1==0");
204     fill_result_histos(zossfpmm, zossfpmmerr, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,datajzb,cut,cuthigh,mcordata,resultmm,"id1==1");
205     }
206    
207 buchmann 1.12 // float ppzossfp, ppzossfperr, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn, ppresult;
208     // if(mcordata==mc) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,cuthigh,mcordata,ppresult);
209     // if(mcordata==data) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(datajzb,jzbpeakerrorData),cut,cuthigh,mcordata,ppresult);
210     // if(mcordata==mcwithsignal) fill_result_histos(ppzossfp, ppzossfperr,ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,cuthigh,mcordata,ppresult);
211 buchmann 1.11
212     float pnzossfp, pnzossfperr, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn, pnresult;
213     if(mcordata==mc) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,cuthigh,mcordata,pnresult);
214     if(mcordata==data) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(datajzb,-jzbpeakerrorData),cut,cuthigh,mcordata,pnresult);
215     if(mcordata==mcwithsignal) fill_result_histos(pnzossfp, pnzossfperr,pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,cuthigh,mcordata,pnresult);
216    
217     float syserr=0;
218     float peakerr=0;
219     float staterr=0;
220     float poissonstaterrup=-999;
221     float poissonstaterrdown=-999;
222 buchmann 1.12 // if(fabs(result-pnresult)>fabs(result-ppresult)) peakerr=fabs(result-pnresult); else peakerr=fabs(result-ppresult);
223 buchmann 1.11
224 buchmann 1.13 float zjetsestimateuncert=0.25;
225 buchmann 1.11 float emuncert=0.5;
226     float emsidebanduncert=0.5;
227     float eemmsidebanduncert=0.5;
228     syserr = (zjetsestimateuncert*zossfn)*(zjetsestimateuncert*zossfn);//first term
229     syserr+= ((zosofp)*(zosofp) + (zosofn)*(zosofn))*(1.0/9)*emuncert*emuncert;//sys err from emu method
230     syserr+= ((sbossfp)*(sbossfp)+(sbossfn)*(sbossfn))*(1.0/9)*eemmsidebanduncert*eemmsidebanduncert; // sys err from eemm sidebands
231     syserr+= ((sbosofp)*(sbosofp)+(sbosofn)*(sbosofn))*(1.0/9)*emsidebanduncert*emsidebanduncert; // sys err from emu sidebands
232     syserr=TMath::Sqrt(syserr);
233    
234     staterr=TMath::Sqrt(zossfn + (1.0/9)*(zosofp+zosofn)+ (1.0/9)*(sbossfp+sbossfn)+ (1.0/9)*(sbosofp+sbosofn));
235     if(dopoisson) advanced_poisson(zossfn,zosofp,zosofn,sbossfp,sbossfn,sbosofp,sbosofn,poissonstaterrdown,poissonstaterrup);
236    
237     float e_to_emu=0.5;
238     float m_to_emu=0.5;
239    
240     if(writeanything) {
241     if(mcordata==mc) dout << " MC :: ";
242     if(mcordata==mcwithsignal) dout << " MC with S :: ";
243     if(mcordata==data) dout << " ";
244     dout << "Observed : " << zossfp << "+-" << zossfperr << endl;
245     if(mcordata==data) dout << " Composition: " << zossfpee << " (ee), " << zossfpmm << " (mm) " << endl;
246    
247     if(mcordata==mc) dout << " MC :: ";
248     if(mcordata==mcwithsignal) dout << " MC with S :: ";
249     if(mcordata==data) dout << " ";
250     dout << "Predicted: " << zossfn << " + (1/3)*(" << zosofp << "-" << zosofn<<") [e&mu]+ (1/3)*(" << sbossfp << "-" << sbossfn<<") [SF,SB]+ (1/3)*(" << sbosofp << "-" << sbosofn<<") [OF,SB] = ";
251     dout << zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn)<< " +/- " << peakerr << " (peak ) " << " +/- " << syserr << " (sys) +/- " << staterr << " (stat)" << endl;
252     if(dopoisson) dout << " stat error with Poisson : +" << poissonstaterrup << " - " << poissonstaterrdown << endl;
253    
254     if(mcordata==data) {
255     dout << " Composition: (ee)" << zossfnee << " + (1/3)*(" << e_to_emu*zosofp << "-" << e_to_emu*zosofn<<") ["<<e_to_emu<<"*e&mu]+ (1/3)*(" << sbossfpee << "-" << sbossfnee<<") [SF,SB]+ (1/3)*(" << sbosofp*e_to_emu << "-" << sbosofn*e_to_emu<<") ["<<e_to_emu<<"*OF,SB] = ";
256     dout << zossfnee + (1.0/3)*(e_to_emu*zosofp-e_to_emu*zosofn)+ (1.0/3)*(sbossfpee-sbossfnee) + (1.0/3)*(sbosofp*e_to_emu - sbosofn*e_to_emu) << endl;
257     dout << " (mm)" << zossfnmm << " + (1/3)*(" << m_to_emu*zosofp << "-" << m_to_emu*zosofn<<") ["<<m_to_emu<<"*e&mu]+ (1/3)*(" << sbossfpmm << "-" << sbossfnmm<<") [SF,SB]+ (1/3)*(" << sbosofp*m_to_emu << "-" << sbosofn*m_to_emu<<") ["<<m_to_emu<<"*OF,SB] = ";
258     dout << zossfnmm + (1.0/3)*(m_to_emu*zosofp - m_to_emu*zosofn) + (1.0/3)*(sbossfpmm -sbossfnmm)+ (1.0/3)*(sbosofp*m_to_emu-sbosofn*e_to_emu) << endl;
259     }
260    
261     if(chatty) {
262     dout << " Pred(ZJets ) \t " << zossfn << endl;
263     dout << " Pred(e&mu;]) \t " << zosofp << "-" << zosofn << " = " << zosofp-zosofn<<endl;
264     dout << " Pred(ossf,sb]) \t " << sbossfp << "-" << sbossfn<<" = "<<sbossfp-sbossfn<<endl;
265     dout << " Pred(osof,sb]) \t " << sbosofp << "-" << sbosofn<<" = "<<sbosofp-sbosofn<<endl;
266     }
267     }//end of writeanything
268    
269     if(mcordata==data) {
270     //store the result!
271     Nobs.push_back(zossfp);
272     Npred.push_back(zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn));
273     float totprederr=0;
274     totprederr+=peakerr*peakerr;
275     totprederr+=syserr*syserr;
276     totprederr+=staterr*staterr;
277     totprederr=TMath::Sqrt(totprederr);
278     Nprederr.push_back(totprederr);
279     }
280     vector<float> resultvector;
281     resultvector.push_back(zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn));
282     resultvector.push_back(TMath::Sqrt(peakerr*peakerr+syserr*syserr+staterr*staterr));
283     resultvector.push_back(zossfp);
284     resultvector.push_back(zossfperr);
285    
286     return resultvector;
287     }
288    
289     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) {
290     get_result_between_two_fixed_jzb_values(cut ,15000, mcjzb,datajzb, mcordata,jzbpeakerrorMC, jzbpeakerrorData, rescan, chatty, dopoisson);
291     }
292 buchmann 1.1
293    
294 buchmann 1.6 void get_result(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbcuts, bool chatty=false, bool dopoisson=false,bool doquick=false) {
295 buchmann 1.1 TCanvas *rescan = new TCanvas("rescan","Result Canvas");
296     for(int icut=0;icut<jzbcuts.size();icut++) {
297 buchmann 1.3 get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,data,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
298 buchmann 1.6 if(!doquick) {
299     get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,mc,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
300     get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,mcwithsignal,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
301     }
302 buchmann 1.1 }
303     delete rescan;
304     }