ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ResultModule.C
Revision: 1.8
Committed: Tue Jul 26 07:47:31 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.7: +41 -20 lines
Log Message: 
Providing numbers in data also separately for ee and mm

File Contents

# Content
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 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, vector<int> sel, samplecollection &sampleC, string addcut="") {
23 string xlabel="JZB [GeV] -- for algoritm internal use only!";
24 bool dosignal=false;
25 if(mcordata==mcwithsignal) {
26 dosignal=true;
27 mcordata=0;
28 }
29 TCut basiccutplus;
30 if(addcut=="") basiccutplus=basiccut;
31 else basiccutplus=basiccut&&addcut.c_str();
32 TH1F *ZOSSFP;
33 TH1F *ZOSOFP;
34 TH1F *ZOSSFN;
35 TH1F *ZOSOFN;
36
37 TH1F *SBOSSFP;
38 TH1F *SBOSOFP;
39 TH1F *SBOSSFN;
40 TH1F *SBOSOFN;
41
42 if(mcordata==mc||mcordata==data||mcordata==mcwithsignal) {
43 ZOSSFP = sampleC.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
44 ZOSOFP = sampleC.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
45 ZOSSFN = sampleC.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
46 ZOSOFN = sampleC.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,dosignal);
47
48 SBOSSFP = sampleC.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
49 SBOSOFP = sampleC.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
50 SBOSSFN = sampleC.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
51 SBOSOFN = sampleC.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,dosignal);
52 } else {
53 //doing signal only!
54 ZOSSFP = sampleC.Draw("ZOSSFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
55 ZOSOFP = sampleC.Draw("ZOSOFP",datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
56 ZOSSFN = sampleC.Draw("ZOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSSF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
57 ZOSOFN = sampleC.Draw("ZOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutmass&&cutOSOF&&cutnJets&&basiccutplus,mcordata,luminosity,sel);
58
59 SBOSSFP = sampleC.Draw("SBOSSFP",datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
60 SBOSOFP = sampleC.Draw("SBOSOFP",datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
61 SBOSSFN = sampleC.Draw("SBOSSFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSSF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
62 SBOSOFN = sampleC.Draw("SBOSOFN","-"+datajzb,1,cut,14000, xlabel, "events",cutOSOF&&cutnJets&&basiccutplus&&sidebandcut,mcordata,luminosity,sel);
63 }
64
65 zossfp=ZOSSFP->Integral();
66 zosofp=ZOSOFP->Integral();
67 zossfn=ZOSSFN->Integral();
68 zosofn=ZOSOFN->Integral();
69
70 sbossfp=SBOSSFP->Integral();
71 sbosofp=SBOSOFP->Integral();
72 sbossfn=SBOSSFN->Integral();
73 sbosofn=SBOSOFN->Integral();
74
75 delete ZOSSFP;
76 delete ZOSOFP;
77 delete ZOSSFN;
78 delete ZOSOFN;
79
80 delete SBOSSFP;
81 delete SBOSOFP;
82 delete SBOSSFN;
83 delete SBOSOFN;
84
85 result = zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn);
86 }
87
88 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, string addcut="") {
89 vector<int> emptyvector;
90 fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,mcordata,result,emptyvector,allsamples,addcut);
91 }
92
93
94 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) {
95 rescan->cd();
96 if(mcordata==data) dout << "Crunching numbers for JZB>" << cut << " : " << endl;
97
98 float zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,result;
99 if(mcordata==mc) fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,mcordata,result);
100 if(mcordata==data) fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,datajzb,cut,mcordata,result);
101 if(mcordata==mcwithsignal) fill_result_histos(zossfp, zosofp, zossfn, zosofn, sbossfp, sbosofp, sbossfn, sbosofn,mcjzb,cut,mcordata,result);
102
103 float zossfpee, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,resultee;
104 float zossfpmm, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,resultmm;
105 if(mcordata==data) {
106 fill_result_histos(zossfpee, zosofpee, zossfnee, zosofnee, sbossfpee, sbosofpee, sbossfnee, sbosofnee,datajzb,cut,mcordata,resultee,"id1==0");
107 fill_result_histos(zossfpmm, zosofpmm, zossfnmm, zosofnmm, sbossfpmm, sbosofpmm, sbossfnmm, sbosofnmm,datajzb,cut,mcordata,resultmm,"id1==1");
108 }
109
110 float ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn, ppresult;
111 if(mcordata==mc) fill_result_histos(ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,mcordata,ppresult);
112 if(mcordata==data) fill_result_histos(ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(datajzb,jzbpeakerrorData),cut,mcordata,ppresult);
113 if(mcordata==mcwithsignal) fill_result_histos(ppzossfp, ppzosofp, ppzossfn, ppzosofn, ppsbossfp, ppsbosofp, ppsbossfn, ppsbosofn,newjzbexpression(mcjzb,jzbpeakerrorMC),cut,mcordata,ppresult);
114
115 float pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn, pnresult;
116 if(mcordata==mc) fill_result_histos(pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,mcordata,pnresult);
117 if(mcordata==data) fill_result_histos(pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(datajzb,-jzbpeakerrorData),cut,mcordata,pnresult);
118 if(mcordata==mcwithsignal) fill_result_histos(pnzossfp, pnzosofp, pnzossfn, pnzosofn, pnsbossfp, pnsbosofp, pnsbossfn, pnsbosofn,newjzbexpression(mcjzb,-jzbpeakerrorMC),cut,mcordata,pnresult);
119
120 float syserr=0;
121 float peakerr=0;
122 float staterr=0;
123 float poissonstaterrup=-999;
124 float poissonstaterrdown=-999;
125 if(fabs(result-pnresult)>fabs(result-ppresult)) peakerr=fabs(result-pnresult); else peakerr=fabs(result-ppresult);
126
127 float zjetsestimateuncert=0.5;
128 float emuncert=0.5;
129 float emsidebanduncert=0.5;
130 float eemmsidebanduncert=0.5;
131 syserr = (zjetsestimateuncert*zossfn)*(zjetsestimateuncert*zossfn);//first term
132 syserr+= ((zosofp)*(zosofp) + (zosofn)*(zosofn))*(1.0/9)*emuncert*emuncert;//sys err from emu method
133 syserr+= ((sbossfp)*(sbossfp)+(sbossfn)*(sbossfn))*(1.0/9)*eemmsidebanduncert*eemmsidebanduncert; // sys err from eemm sidebands
134 syserr+= ((sbosofp)*(sbosofp)+(sbosofn)*(sbosofn))*(1.0/9)*emsidebanduncert*emsidebanduncert; // sys err from emu sidebands
135 syserr=TMath::Sqrt(syserr);
136
137 staterr=TMath::Sqrt(zossfn + (1.0/9)*(zosofp+zosofn)+ (1.0/9)*(sbossfp+sbossfn)+ (1.0/9)*(sbosofp+sbosofn));
138 if(dopoisson) advanced_poisson(zossfn,zosofp,zosofn,sbossfp,sbossfn,sbosofp,sbosofn,poissonstaterrdown,poissonstaterrup);
139
140 float e_to_emu=0.5;
141 float m_to_emu=0.5;
142
143 if(mcordata==mc) dout << " MC :: ";
144 if(mcordata==mcwithsignal) dout << " MC with S :: ";
145 if(mcordata==data) dout << " ";
146 dout << "Observed : " << zossfp << endl;
147 if(mcordata==data) dout << " Composition: " << zossfpee << " (ee), " << zossfpmm << " (mm) " << endl;
148
149 if(mcordata==mc) dout << " MC :: ";
150 if(mcordata==mcwithsignal) dout << " MC with S :: ";
151 if(mcordata==data) dout << " ";
152 dout << "Predicted: " << zossfn << " + (1/3)*(" << zosofp << "-" << zosofn<<") [e&mu]+ (1/3)*(" << sbossfp << "-" << sbossfn<<") [SF,SB]+ (1/3)*(" << sbosofp << "-" << sbosofn<<") [OF,SB] = ";
153 dout << zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn)<< " +/- " << peakerr << " (peak ) " << " +/- " << syserr << " (sys) +/- " << staterr << " (stat)" << endl;
154 if(dopoisson) dout << " stat error with Poisson : +" << poissonstaterrup << " - " << poissonstaterrdown << endl;
155
156 if(mcordata==data) {
157 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] = ";
158 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;
159 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] = ";
160 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;
161 }
162
163 if(chatty) {
164 dout << " Pred(ZJets ) \t " << zossfn << endl;
165 dout << " Pred(e&mu;]) \t " << zosofp << "-" << zosofn << " = " << zosofp-zosofn<<endl;
166 dout << " Pred(ossf,sb]) \t " << sbossfp << "-" << sbossfn<<" = "<<sbossfp-sbossfn<<endl;
167 dout << " Pred(osof,sb]) \t " << sbosofp << "-" << sbosofn<<" = "<<sbosofp-sbosofn<<endl;
168 }
169
170 if(mcordata==data) {
171 //store the result!
172 Nobs.push_back(zossfp);
173 Npred.push_back(zossfn + (1.0/3)*(zosofp-zosofn)+ (1.0/3)*(sbossfp-sbossfn)+ (1.0/3)*(sbosofp-sbosofn));
174 float totprederr=0;
175 totprederr+=peakerr*peakerr;
176 totprederr+=syserr*syserr;
177 totprederr+=staterr*staterr;
178 totprederr=TMath::Sqrt(totprederr);
179 Nprederr.push_back(totprederr);
180 }
181 }
182
183
184 void get_result(string mcjzb, string datajzb, float jzbpeakerrordata, float jzbpeakerrormc, vector<float> jzbcuts, bool chatty=false, bool dopoisson=false,bool doquick=false) {
185 TCanvas *rescan = new TCanvas("rescan","Result Canvas");
186 for(int icut=0;icut<jzbcuts.size();icut++) {
187 get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,data,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
188 if(!doquick) {
189 get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,mc,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
190 get_result_above_one_fixed_jzb_value(jzbcuts[icut],mcjzb,datajzb,mcwithsignal,jzbpeakerrormc,jzbpeakerrordata,rescan,chatty,dopoisson);
191 }
192 }
193 delete rescan;
194 }