ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Plotting/Modules/ResultModule.C
Revision: 1.6
Committed: Tue Jul 19 15:01:56 2011 UTC (13 years, 9 months ago) by buchmann
Content type: text/plain
Branch: MAIN
Changes since 1.5: +59 -15 lines
Log Message: 
Adapted the result module to be conform with limit treatment and permit calculation of expected and observed yields for signal

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