ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/Development/Plotting/Modules/DustyModule.C
Revision: 1.2
Committed: Wed Sep 5 20:17:20 2012 UTC (12 years, 8 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +195 -0 lines
Log Message: 
Added experimental met z prediction to dusty module; the refined version will be uploaded to MetPlotting in a couple of seconds

File Contents

# User Rev Content
1 buchmann 1.1 /****
2    
3     Off peak status (RestrictToMassPeak) :
4    
5     THIS MODULE WILL NOT BE ADAPTED AS IT ONLY CONTAINS OUTDATED FUNCTIONS!
6    
7    
8     . Necessary adaptations identified
9     . Started working on necessary adaptations
10     . Necessary adaptations implemented
11     . Necessary adaptations tested
12    
13    
14    
15     ****/
16    
17     /*
18    
19     Contains functions that were used in the past but are out now ...
20    
21     */
22     #include <iostream>
23     #include <vector>
24     #include <sys/stat.h>
25    
26     #include <TCut.h>
27     #include <TROOT.h>
28     #include <TCanvas.h>
29     #include <TMath.h>
30     #include <TColor.h>
31     #include <TPaveText.h>
32     #include <TRandom.h>
33     #include <TH1.h>
34     #include <TH2.h>
35     #include <TF1.h>
36     #include <TSQLResult.h>
37     #include <TProfile.h>
38    
39     string centralupdownname(int num) {
40     if (num==0) return "central";
41     if (num==1) return "down";
42     if (num==2) return "up";
43     return "centralupdownnameERROR";
44     }
45    
46     string dataormc(int isdata) {
47     if(isdata) return "Data";
48     else return "MC";
49     }
50    
51     int histocounter=0;
52     string give_histo_number(string basename="", int isdata=-1, int centralcounter=-1) {
53     histocounter++;
54     stringstream histocounternum;
55     if(isdata==-1) {
56     if(basename=="") histocounternum << histocounter;
57     else histocounternum << basename << "_" << histocounter;
58     }
59     else {
60     histocounternum << basename << "_" << dataormc(isdata) << "_" << centralupdownname(centralcounter);
61     }
62     return histocounternum.str();
63     }
64    
65     int central=0;
66     int down=1;
67     int up=2;
68    
69     void crunch_the_numbers(TH1F *RcorrJZBeemmop[3][2],TH1F *RcorrJZBeeop[3][2],TH1F *RcorrJZBmmop[3][2],TH1F *LcorrJZBeemmop[3][2],TH1F *LcorrJZBeeop[3][2],TH1F *LcorrJZBmmop[3][2],TH1F *RcorrJZBemop[3][2],TH1F *LcorrJZBemop[3][2],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],float predicted[3][2][20],float observed[3][2][20],int isdata,vector<float> jzbcuts,float eemmmumu[2][3][20]) {
70     // dout << "Crunching the numbers for is_data=" << isdata << endl;
71     for(int k=0;k<20;k++) {for(int i=0;i<2;i++){for(int j=0;j<3;j++) {eemmmumu[i][j][k]=0;}}}
72     for(int icut=0;icut<jzbcuts.size()-1;icut++) {//do calculation for each JZB cut
73     for(int ibin=1;ibin<jzbcuts.size();ibin++) {//and to do that, we need to consider each bin.
74     if(RcorrJZBeemmop[0][isdata]->GetBinLowEdge(ibin)<jzbcuts[icut]) continue;
75     // if(icut==0) dout << "Predicted+=" << LcorrJZBeemmop[central][isdata]->GetBinContent(ibin) << " eemm,L" << endl;
76     predicted[central][isdata][icut]+=LcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
77     predicted[up][isdata][icut]+=LcorrJZBeemmop[up][isdata]->GetBinContent(ibin);
78     predicted[down][isdata][icut]+=LcorrJZBeemmop[down][isdata]->GetBinContent(ibin);
79    
80     // if(icut==0) dout << "Predicted+=" << RcorrJZBemop[central][isdata]->GetBinContent(ibin) << " em, R" << endl;
81     predicted[central][isdata][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
82     predicted[up][isdata][icut]+=RcorrJZBemop[up][isdata]->GetBinContent(ibin);
83     predicted[down][isdata][icut]+=RcorrJZBemop[down][isdata]->GetBinContent(ibin);
84    
85     // if(icut==0) dout << "Predicted-=" << LcorrJZBemop[central][isdata]->GetBinContent(ibin) <<" em, L" << endl;
86     predicted[central][isdata][icut]-=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
87     predicted[up][isdata][icut]-=LcorrJZBemop[up][isdata]->GetBinContent(ibin);
88     predicted[down][isdata][icut]-=LcorrJZBemop[down][isdata]->GetBinContent(ibin);
89    
90     // if(icut==0) dout << "Observed=" << RcorrJZBeemmop[central][isdata]->GetBinContent(ibin) << "eemm, R" << endl;
91     observed[central][isdata][icut]+=RcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
92     observed[up][isdata][icut]+=RcorrJZBeemmop[up][isdata]->GetBinContent(ibin);
93     observed[down][isdata][icut]+=RcorrJZBeemmop[down][isdata]->GetBinContent(ibin);
94    
95     zjetsestimate[central][isdata][icut]+=LcorrJZBeemmop[central][isdata]->GetBinContent(ibin);
96     ttbarestimate[central][isdata][0][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
97     ttbarestimate[central][isdata][1][icut]+=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
98    
99     eemmmumu[0][0][icut]+=LcorrJZBeeop[central][isdata]->GetBinContent(ibin);
100     eemmmumu[0][1][icut]+=LcorrJZBmmop[central][isdata]->GetBinContent(ibin);
101     eemmmumu[0][2][icut]+=LcorrJZBemop[central][isdata]->GetBinContent(ibin);
102    
103     eemmmumu[1][0][icut]+=RcorrJZBeeop[central][isdata]->GetBinContent(ibin);
104     eemmmumu[1][1][icut]+=RcorrJZBmmop[central][isdata]->GetBinContent(ibin);
105     eemmmumu[1][2][icut]+=RcorrJZBemop[central][isdata]->GetBinContent(ibin);
106    
107     }//end of ibin loop
108     }//end of icut loop
109     }
110    
111     void make_table(float eemmmumu[2][3][20],int icut,float JZBcut) {
112     vector< vector <string> > entries;
113     vector <string> line;
114     line.push_back("");
115     line.push_back("eemm (ee/mm)");
116     line.push_back("em");
117     entries.push_back(line);
118     line.clear();
119     line.push_back("JZB>"+any2string(JZBcut));
120     line.push_back(any2string(eemmmumu[1][0][icut]+eemmmumu[1][1][icut])+"("+any2string(eemmmumu[1][0][icut])+"/"+any2string(eemmmumu[1][1][icut])+")");
121     line.push_back(any2string(eemmmumu[1][2][icut]));
122     entries.push_back(line);
123     line.clear();
124     line.push_back("JZB<-"+any2string(JZBcut));
125     line.push_back(any2string(eemmmumu[0][0][icut]+eemmmumu[0][1][icut])+"("+any2string(eemmmumu[0][0][icut])+"/"+any2string(eemmmumu[0][1][icut])+")");
126     line.push_back(any2string(eemmmumu[0][2][icut]));
127     entries.push_back(line);
128     make_nice_jzb_table(entries);
129     /*
130     dout << " \t\t | \t eemm (ee/mm) \t | \t em" << endl;
131     dout << "JZB>peak\t | \t " << eemmmumu[1][0]+eemmmumu[1][1] << "(" << eemmmumu[1][0] << "/" << eemmmumu[1][1] << ")\t | \t" << eemmmumu[1][2] << endl;
132     dout << "JZB<peak\t | \t " << eemmmumu[0][0]+eemmmumu[0][1] << "(" << eemmmumu[0][0] << "/" << eemmmumu[0][1] << ")\t | \t" << eemmmumu[0][2] << endl;
133     dout << endl;
134     */
135     }
136    
137     void present_result(vector<float> &jzbcuts,float predicted[3][2][20],float observed[3][2][20],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],int icut, float eemmmumu[2][3][20]) {
138    
139     //blublu
140     dout << endl << endl;
141     dout << "--------------------------------------------------------------------------------------" << endl;
142     dout << "DATA: " << endl;
143     dout << " For JZB>" << jzbcuts[icut];
144     float max,downvar,upvar;
145     stringstream printtitle;
146     printtitle << "JZB>" << jzbcuts[icut] << " (data)";
147     ComputePoissonError(zjetsestimate[central][data][icut],ttbarestimate[central][data][0][icut],ttbarestimate[central][data][1][icut],max,downvar,upvar,printtitle.str());
148     //dout << " Predicted: " << print_range(predicted[central][data][icut],predicted[up][data][icut],predicted[down][data][icut]) << " (stat) +" << 0.2*zjetsestimate[central][data][icut] << "-" << 0.4*zjetsestimate[central][data][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
149     float pred=predicted[central][data][icut];
150     float sysP=abs(predicted[up][data][icut]-pred);
151     float sysN=abs(pred-predicted[down][data][icut]);
152     sysN = sysN + pred*(1-fitresultconstdata);//fitresultconst comes from the fit in the 0-30 GeV range!
153     dout << " Predicted: " << pred << "+" << statErrorP(pred) << "-" << statErrorN(pred) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
154     //0.2*zjetsestimate[central][data][icut] << "-" << 0.4*zjetsestimate[central][data][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
155     dout << " Details: ZJetsEstimate= " << zjetsestimate[central][data][icut] << ", TTbar estimate=" << ttbarestimate[central][data][0][icut]-ttbarestimate[central][data][1][icut] << " [" << ttbarestimate[central][data][0][icut] << " R -" << ttbarestimate[central][data][1][icut] << " L ]" << endl;
156     dout << " For JZB>" << jzbcuts[icut] << " Observed: " << observed[central][data][icut] << endl;
157     dout << "TABLE:" << endl;
158     make_table(eemmmumu,icut,jzbcuts[icut]);
159     dout << "MC: " << endl;
160     printtitle.str("");
161     printtitle<<"JZB>"<<jzbcuts[icut]<<" (MC)";
162     ComputePoissonError(zjetsestimate[central][mc][icut],ttbarestimate[central][mc][0][icut],ttbarestimate[central][mc][1][icut],max,downvar,upvar,printtitle.str());
163     pred=predicted[central][mc][icut];
164     sysP=abs(predicted[up][mc][icut]-pred);
165     sysN=abs(pred-predicted[down][mc][icut]);
166     sysN = sysN + pred*(1-fitresultconstmc);//fitresultconst comes from the fit in the 0-30 GeV range!
167     dout << " Predicted: " << pred << "+" << statErrorP(pred) << "-" << statErrorN(pred) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
168     // dout << " Predicted: " << print_range(pred,statErrorP(pred),statErrorN(pred)) << " (stat) +" << sysP << " - " << sysN << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
169     //dout << " Predicted: " << print_range(predicted[central][mc][icut],predicted[down][mc][icut],predicted[up][mc][icut]) << " (stat) +" << 0.2*zjetsestimate[central][mc][icut] << "-" << 0.4*zjetsestimate[central][mc][icut] << " (sys) " << " using Poisson: " << print_range(max,max+upvar,max-downvar) << ")" << endl;
170     //print_range(predicted[central][mc][icut],predicted[central][mc][icut]+upvar,predicted[central][mc][icut]-downvar) << ")" << endl;
171     dout << " Details: ZJetsEstimate= " << zjetsestimate[central][mc][icut] << ", TTbar estimate=" << ttbarestimate[central][mc][0][icut]-ttbarestimate[central][mc][1][icut] << " [" << ttbarestimate[central][mc][0][icut] << " R -" << ttbarestimate[central][mc][1][icut] << " L ] )" << endl;
172     dout << " For JZB>" << jzbcuts[icut] << " Observed: " << observed[central][mc][icut] << endl;
173     dout << endl << endl;
174     }
175    
176     void calculate_predicted_and_observed_eemm(float MCPeak,float MCPeakError,float DataPeak,float DataPeakError,vector<float> jzbcuts) { // DO NOT CHANGE THIS TO ... &jzbcuts !!! we add one!
177     jzbcuts.push_back(14000);
178    
179     /*
180     * We want the following numbers: For all JZB cuts, we want: ee,mm,em; observed, predicted. we want final results with errors from the peak and statistical error.
181     * How to accomplish this; we draw histos for MC&data, once with the peak (to extract obs/pred), once with peak+sigma, once with peak-sigma. this gives us the peak error.
182     * Statistical error: Comes from Poissonian errors ... for this we have a special little macro that will give us the value for three given parameters (cool, right?)
183     */
184    
185     string jzbpeak[3][2];
186     jzbpeak[central][data]=give_jzb_expression(DataPeak,data);
187     jzbpeak[up][data]=give_jzb_expression(DataPeak+DataPeakError,data);
188     jzbpeak[down][data]=give_jzb_expression(DataPeak-DataPeakError,data);
189     jzbpeak[central][mc]=give_jzb_expression(MCPeak,mc);
190     jzbpeak[up][mc]=give_jzb_expression(MCPeak+MCPeakError,mc);
191     jzbpeak[down][mc]=give_jzb_expression(MCPeak-MCPeakError,mc);
192    
193     TH1F *RcorrJZBeemmop[3][2];
194     TH1F *LcorrJZBeemmop[3][2];
195     TH1F *RcorrJZBeeop[3][2];
196     TH1F *LcorrJZBeeop[3][2];
197     TH1F *RcorrJZBmmop[3][2];
198     TH1F *LcorrJZBmmop[3][2];
199     TH1F *RcorrJZBemop[3][2];
200     TH1F *LcorrJZBemop[3][2];
201     float zjetsestimate[3][2][20];
202     float predicted[3][2][20];
203     float observed[3][2][20];
204     float ttbarestimate[3][2][2][20];
205     float eemmmumu[2][3][20];
206    
207    
208     for(int isdata=0;isdata<=1;isdata++) {
209     for(int centraldownup=0;centraldownup<=2;centraldownup++) {
210     RcorrJZBeemmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBeemmop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,isdata, luminosity);
211     LcorrJZBeemmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBeemmop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,isdata, luminosity);
212    
213     RcorrJZBemop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBemop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,isdata, luminosity);
214     LcorrJZBemop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBemop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,isdata,luminosity);
215     RcorrJZBemop[centraldownup][isdata]->Add(allsamples.Draw(give_histo_number("RcorrJZBemopSB",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,isdata, luminosity));
216     LcorrJZBemop[centraldownup][isdata]->Add(allsamples.Draw(give_histo_number("LcorrJZBemopSB",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", sidebandcut&&cutOSOF&&cutnJets,isdata,luminosity));
217     RcorrJZBemop[centraldownup][isdata]->Add(allsamples.Draw(give_histo_number("RcorrJZBeemmopSB",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,isdata, luminosity));
218     LcorrJZBemop[centraldownup][isdata]->Add(allsamples.Draw(give_histo_number("LcorrJZBeemmopSB",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", sidebandcut&&cutOSSF&&cutnJets,isdata,luminosity));
219     LcorrJZBemop[centraldownup][isdata]->Scale(1.0/3);
220     RcorrJZBemop[centraldownup][isdata]->Scale(1.0/3);
221     RcorrJZBeeop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBeeop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"id1==0",isdata, luminosity);
222     LcorrJZBeeop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBeeop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"id1==0",isdata, luminosity);
223     RcorrJZBmmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("RcorrJZBmmop",isdata,centraldownup),jzbpeak[centraldownup][isdata].c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"id1==1",isdata, luminosity);
224     LcorrJZBmmop[centraldownup][isdata] = allsamples.Draw(give_histo_number("LcorrJZBmmop",isdata,centraldownup),("-"+jzbpeak[centraldownup][isdata]).c_str(),jzbcuts, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"id1==1",isdata, luminosity);
225     for(int icut=0;icut<=jzbcuts.size();icut++) {
226     zjetsestimate[centraldownup][isdata][icut]=0;
227     predicted[centraldownup][isdata][icut]=0;
228     observed[centraldownup][isdata][icut]=0;
229     ttbarestimate[centraldownup][isdata][0][icut]=0;
230     ttbarestimate[centraldownup][isdata][1][icut]=0;
231     }
232     }//end of for central,up,down loop
233     crunch_the_numbers(RcorrJZBeemmop,RcorrJZBeeop,RcorrJZBmmop,LcorrJZBeemmop,LcorrJZBeeop,LcorrJZBmmop,RcorrJZBemop,LcorrJZBemop,zjetsestimate,ttbarestimate,predicted,observed,isdata,jzbcuts,eemmmumu);
234    
235     // void crunch_the_numbers(TH1F *RcorrJZBeemmop[3][2],TH1F *RcorrJZBeeop[3][2],TH1F *RcorrJZBmmop[3][2],TH1F *LcorrJZBeemmop[3][2],TH1F *LcorrJZBeeop[3][2],TH1F *LcorrJZBmmop[3][2],TH1F *RcorrJZBemop[3][2],TH1F *LcorrJZBemop[3][2],float zjetsestimate[3][2][20],float ttbarestimate[3][2][2][20],float predicted[3][2][20],float observed[3][2][20],int isdata,vector<float> jzbcuts,float eemmmumu[2][3]) {
236    
237     }//end of is data loop
238    
239    
240     dout << "We obtain the following results: " << endl;
241     for(int icut=0;icut<jzbcuts.size()-1;icut++) {
242     present_result(jzbcuts,predicted,observed,zjetsestimate,ttbarestimate,icut,eemmmumu);
243     }
244    
245    
246     }
247    
248    
249    
250    
251    
252    
253    
254    
255    
256    
257    
258    
259    
260     void alpha_scan_susy_space(string mcjzb, string datajzb) {
261     TCanvas *c3 = new TCanvas("c3","c3");
262     vector<float> binning;
263     binning=allsamples.get_optimal_binsize(mcjzb,cutmass&&cutOSSF&&cutnJets,20,50,800);
264     float arrbinning[binning.size()];
265     for(int i=0;i<binning.size();i++) arrbinning[i]=binning[i];
266     TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
267     puredata->SetMarkerSize(DataMarkerSize);
268     TH1F *allbgs = allsamples.Draw("allbgs", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
269     TH1F *allbgsb = allsamples.Draw("allbgsb", "-"+datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
270     TH1F *allbgsc = allsamples.Draw("allbgsc", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,luminosity);
271     allbgs->Add(allbgsb,-1);
272     allbgs->Add(allbgsc);
273     int ndata=puredata->Integral();
274     ofstream myfile;
275     myfile.open ("susyscan_log.txt");
276     TFile *susyscanfile = new TFile("/scratch/fronga/SMS/T5z_SqSqToQZQZ_38xFall10.root");
277     TTree *suevents = (TTree*)susyscanfile->Get("events");
278     TH2F *exclusionmap = new TH2F("exclusionmap","",20,0,500,20,0,1000);
279     TH2F *exclusionmap1s = new TH2F("exclusionmap1s","",20,0,500,20,0,1000);
280     TH2F *exclusionmap2s = new TH2F("exclusionmap2s","",20,0,500,20,0,1000);
281     TH2F *exclusionmap3s = new TH2F("exclusionmap3s","",20,0,500,20,0,1000);
282    
283     susy_scan_axis_labeling(exclusionmap);
284     susy_scan_axis_labeling(exclusionmap1s);
285     susy_scan_axis_labeling(exclusionmap2s);
286     susy_scan_axis_labeling(exclusionmap3s);
287    
288     Int_t MyPalette[100];
289     Double_t r[] = {0., 0.0, 1.0, 1.0, 1.0};
290     Double_t g[] = {0., 0.0, 0.0, 1.0, 1.0};
291     Double_t b[] = {0., 1.0, 0.0, 0.0, 1.0};
292     Double_t stop[] = {0., .25, .50, .75, 1.0};
293     Int_t FI = TColor::CreateGradientColorTable(5, stop, r, g, b, 100);
294     for (int i=0;i<100;i++) MyPalette[i] = FI+i;
295    
296     gStyle->SetPalette(100, MyPalette);
297    
298     for(int m23=50;m23<500;m23+=25) {
299     for (int m0=(2*(m23-50)+150);m0<=1000;m0+=50)
300     {
301     c3->cd();
302     stringstream drawcondition;
303     drawcondition << "pfJetGoodNum>=3&&(TMath::Abs(masses[0]-"<<m0<<")<10&&TMath::Abs(masses[2]-masses[3]-"<<m23<<")<10)&&mll>5&&id1==id2";
304     TH1F *puresignal = new TH1F("puresignal","puresignal",binning.size()-1,arrbinning);
305     TH1F *puresignall= new TH1F("puresignall","puresignal",binning.size()-1,arrbinning);
306     stringstream drawvar,drawvar2;
307     drawvar<<mcjzb<<">>puresignal";
308     drawvar2<<"-"<<mcjzb<<">>puresignall";
309     suevents->Draw(drawvar.str().c_str(),drawcondition.str().c_str());
310     suevents->Draw(drawvar2.str().c_str(),drawcondition.str().c_str());
311     if(puresignal->Integral()<60) {
312     delete puresignal;
313     continue;
314     }
315     puresignal->Add(puresignall,-1);//we need to correct for the signal contamination - we effectively only see (JZB>0)-(JZB<0) !!
316     puresignal->Scale(ndata/(20*puresignal->Integral()));//normalizing it to 5% of the data
317     stringstream saveas;
318     saveas<<"Model_Scan/m0_"<<m0<<"__m23_"<<m23;
319     dout << "PLEASE KEEP IN MIND THAT SIGNAL CONTAMINATION IS NOT REALLY TAKEN CARE OF YET DUE TO LOW STATISTICS! SHOULD BE SOMETHING LIKE THIS : "<< endl;
320     // TH1F *signalpredlo = allsamples.Draw("signalpredlo", "-"+mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
321     // TH1F *signalpredro = allsamples.Draw("signalpredro", mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
322     // TH1F *puredata = allsamples.Draw("puredata", datajzb,binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,luminosity);
323     // signalpred->Add(signalpredlo,-1);
324     // signalpred->Add(signalpredro);
325     // puresignal->Add(signalpred,-1);//subtracting signal contamination
326     //---------------------
327     // dout << "(m0,m23)=("<<m0<<","<<m23<<") contains " << puresignal->Integral() << endl;
328     // TH1F *puresignal = allsamples.Draw("puresignal",mcjzb, binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,mc, luminosity,allsamples.FindSample("LM4"));
329     vector<float> results=establish_upper_limits(puredata,allbgs,puresignal,saveas.str(),myfile);
330     if(results.size()==0) {
331     delete puresignal;
332     continue;
333     }
334     exclusionmap->Fill(m23,m0,results[0]);
335     exclusionmap1s->Fill(m23,m0,results[1]);
336     exclusionmap2s->Fill(m23,m0,results[2]);
337     exclusionmap3s->Fill(m23,m0,results[3]);
338     delete puresignal;
339     dout << "(m0,m23)=("<<m0<<","<<m23<<") : 3 sigma at " << results[3] << endl;
340     }
341     }//end of model scan for loop
342    
343     dout << "Exclusion Map contains" << exclusionmap->Integral() << " (integral) and entries: " << exclusionmap->GetEntries() << endl;
344     c3->cd();
345     exclusionmap->Draw("CONTZ");
346     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_Mean_values");
347     exclusionmap->Draw("COLZ");
348     CompleteSave(c3,"Model_Scan/COL/Model_Scan_Mean_values");
349    
350     exclusionmap1s->Draw("CONTZ");
351     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_1sigma_values");
352     exclusionmap1s->Draw("COLZ");
353     CompleteSave(c3,"Model_Scan/COL/Model_Scan_1sigma_values");
354    
355     exclusionmap2s->Draw("CONTZ");
356     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_2sigma_values");
357     exclusionmap2s->Draw("COLZ");
358     CompleteSave(c3,"Model_Scan/COL/Model_Scan_2sigma_values");
359    
360     exclusionmap3s->Draw("CONTZ");
361     CompleteSave(c3,"Model_Scan/CONT/Model_Scan_3sigma_values");
362     exclusionmap3s->Draw("COLZ");
363     CompleteSave(c3,"Model_Scan/COL/Model_Scan_3sigma_values");
364    
365     TFile *exclusion_limits = new TFile("exclusion_limits.root","RECREATE");
366     exclusionmap->Write();
367     exclusionmap1s->Write();
368     exclusionmap2s->Write();
369     exclusionmap3s->Write();
370     exclusion_limits->Close();
371     susyscanfile->Close();
372    
373     myfile.close();
374     }
375    
376    
377    
378    
379     //-------------------------------------------------------------
380    
381     /*
382    
383    
384     not in use anymore:
385    
386     -------------------------------------<Spring11 - used up to Jul 14>----------------------------------------------
387    
388     // allsamples.AddSample("/scratch/buchmann/AllData_Jun10___DCS_648pb_MoreTriggers4_DCS.root","Data",1,1,true,false,0,kBlack);
389     // allsamples.AddSample("/scratch/buchmann/AllData_Jun10___486pb_MoreTriggers3.root","Data",1,1,true,false,0,kBlack);
390     //allsamples.AddSample("/scratch/buchmann/AllData_Jun17_160404-166861_real_complete_2.root","Data",1,1,true,false,0,kBlack);
391     allsamples.AddSample("/scratch/buchmann/ntuples2011/data/AllData_Jul6_1078pb.root","Data",1,1,true,false,0,kBlack);
392    
393     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/TTJets_TuneZ2_7TeV-madgraph-tauola.root","t#bar{t}",totEventsTTbar,TTbarCrossSection,false,false,1,ttbar_color);//spring11
394     allsamples.AddSample("/scratch/fronga/MC_v1.29/TTJets_TuneZ2_7TeV-madgraph-tauola_Summer11_3.root","TTbar",totEventsTTbar,TTbarCrossSection,false,false,1,ttbar_color);//summer11
395     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets",totEventsWJets,WJetsCrossSection,false,false,3,wjets_color); //spring11
396     allsamples.AddSample("/scratch/fronga/MC_v1.29/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets",totEventsWJets,WJetsCrossSection,false,false,3,wjets_color);//summer
397     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root","DiBosons",totEventsVVJets,VVJetsCrossSection,false,false,4,diboson_color);
398     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/TToBLNu_TuneZ2_s-channel_7TeV-madgraph_2.root","SingleTop",totEventsSingleTopS,SingleTopSCrossSection,false,false,2,singletop_color);
399     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopT,SingleTopTCrossSection,false,false,2,singletop_color);
400     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopU,SingleTopUCrossSection,false,false,2,singletop_color);
401     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/ZinvisibleJets_7TeV-madgraph.root","Z nunu",totEventsZnunu,ZnunuCrossSection,false,false,4,diboson_color);
402     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/QCD_TuneD6T_HT-100To250_7TeV-madgraph.root","QCD",totEventsQCD100to250,QCD100to250CrossSection,false,false,5,qcd_color);
403     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/QCD_TuneD6T_HT-250To500_7TeV-madgraph.root","QCD",totEventsQCD250to500,QCD250to500CrossSection,false,false,5,qcd_color);
404     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/QCD_TuneD6T_HT-500To1000_7TeV-madgraph.root","QCD",totEventsQCD500to1000,QCD500to1000CrossSection,false,false,5,qcd_color);
405     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/QCD_TuneD6T_HT-1000_7TeV-madgraph.root","QCD",totEventsQCD1000toInf,QCD1000toInfCrossSection,false,false,5,qcd_color);
406     // allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root","ZJets",totEventsZjets,ZJetsCrossSection,false,false,6,dy_color);//spring11
407     allsamples.AddSample("/scratch/fronga/MC_v1.29/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola_Summer11.root","Z+Jets",totEventsZjets,ZJetsCrossSection,false,false,6,dy_color);//summer11
408     allsamples.AddSample("/scratch/buchmann/MCSpring2011PU__inclPFvarss__PURW_160404_166861/LM4_SUSY_sftsht_7TeV-pythia6.root","LM4",218536,1.879,false,true,7,lm_color);
409    
410     -------------------------------------</Spring11 - used up to Jul 14>----------------------------------------------
411    
412    
413    
414    
415    
416     //how this works: AddSample("FilePath,"Sample name",Nevents,XS,is_data,is_signal,groupindex
417     allsamples.AddSample("/scratch/buchmann/AllDataNov4.root","Data",1,1,true,false,0,kBlack);
418     allsamples.AddSample("/scratch/buchmann/CutAnalysis/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets to LNu",14350756,31314,false,false,2,kGreen);
419     // allsamples.AddSample("/scratch/buchmann/CutAnalysis/DYJetsToLL_TuneD6T_M-10To50_7TeV-madgraph-tauola.root","DY,M10To50",186891,400,false,false,3);
420     allsamples.AddSample("/scratch/buchmann/CutAnalysis/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root","DY,M50",2552567,3048,false,false,3,kYellow);
421     allsamples.AddSample("/scratch/buchmann/CutAnalysis/TTJets_TuneZ2_7TeV-madgraph-tauola.root","TTJets",1144028,157.5,false,false,4,kMagenta-1);
422     allsamples.AddSample("/scratch/buchmann/CutAnalysis/TToBLNu_TuneZ2_s-channel_7TeV-madgraph.root","T to blnu,s",489472,1.49,false,false,5,kBlue);
423     allsamples.AddSample("/scratch/buchmann/CutAnalysis/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root","T to blnu,t",477610,20.54,false,false,5,kBlue);
424     allsamples.AddSample("/scratch/buchmann/CutAnalysis/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root","T to blnu,u",477599,10.6,false,false,5,kBlue);
425     allsamples.AddSample("/scratch/buchmann/CutAnalysis/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root","VVJets to 4L",509072,4.8,false,false,1,kGreen);
426     allsamples.AddSample("/scratch/buchmann/CutAnalysis/LM4_SUSY_sftsht_7TeV-pythia6.root","20xLM4",218536,1.879*100,false,true,6,kMagenta);
427    
428    
429    
430     /*
431     Int_t fermiwjets = TColor::GetColor("#75f075");
432     Int_t ferminonW = TColor::GetColor("#7d99d1");
433     Int_t fermiwwwz = TColor::GetColor("#2b26ac");
434     Int_t fermizjets = TColor::GetColor("#0044ff");
435     Int_t fermitop = TColor::GetColor("#353535");
436     Int_t fermihiggs = TColor::GetColor("#b30000");
437    
438    
439     dy_color = fermiwjets;
440     wjets_color = ferminonW;
441     ttbar_color = fermizjets;
442     singletop_color = fermitop;
443     qcd_color = nice_orange;
444     diboson_color = fermiwwwz;
445     lm_color = fermihiggs;
446     */
447    
448    
449     /*
450     Int_t dy_color = TColor::GetColor("#ff1300");
451     Int_t wjets_color = TColor::GetColor("#ff7a00");
452     Int_t ttbar_color = TColor::GetColor("#03899c");
453     Int_t singletop_color = TColor::GetColor("#00c322");
454     Int_t qcd_color = TColor::GetColor("#4a11ae");
455     Int_t diboson_color = TColor::GetColor("#8f04a8");
456     Int_t lm_color = TColor::GetColor("#dcf900");
457    
458     //fugly set of colors:
459     dy_color = TColor::GetColor("#ff1300");
460     wjets_color = TColor::GetColor("#ff7a00");
461     ttbar_color = TColor::GetColor("#03899c");
462     singletop_color = TColor::GetColor("#00c322");
463     qcd_color = TColor::GetColor("#4a11ae");
464     diboson_color = TColor::GetColor("#8f04a8");
465     lm_color = TColor::GetColor("#dcf900");
466    
467     //twilight
468     dy_color = TColor::GetColor("#141414");
469     wjets_color = TColor::GetColor("#f8f8f8");
470     ttbar_color = TColor::GetColor("#cf6a4c");
471     singletop_color = TColor::GetColor("#cda869");
472     qcd_color = TColor::GetColor("#e9c062");
473     diboson_color = TColor::GetColor("#5f5a60");
474     lm_color = TColor::GetColor("#7587a6");
475     //http://www.colorschemer.com/schemes/viewscheme.php?id=9126
476    
477     */
478    
479     /* //how this works: AddSample("FilePath,"Sample name",Nevents,XS,is_data,is_signal,groupindex
480     allsamples.AddSample("/scratch/fronga/JSON_160404-163869_v1.29/AllCertified191_v1.29.root","Data",1,1,true,false,0,kBlack);
481     allsamples.AddSample("/scratch/fronga/MC_v1.29/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root","ZJets",totEventsZjets,ZJetsCrossSection,false,false,3,kYellow);
482     allsamples.AddSample("/scratch/fronga/MC_v1.29/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets",totEventsWJets,WJetsCrossSection,false,false,2,kGray);
483     allsamples.AddSample("/scratch/fronga/MC_v1.29/TTJets_TuneZ2_7TeV-madgraph-tauola.root","t#bar{t}",totEventsTTbar,TTbarCrossSection,false,false,4,kMagenta-1);
484     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_s-channel_7TeV-madgraph_2.root","SingleTop",totEventsSingleTopS,SingleTopSCrossSection,false,false,5,kBlue);
485     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopT,SingleTopTCrossSection,false,false,5,kBlue);
486     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopU,SingleTopUCrossSection,false,false,5,kBlue);
487     allsamples.AddSample("/scratch/fronga/MC_v1.29/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root","Dibosons",totEventsVVJets,VVJetsCrossSection,false,false,1,kGreen);
488     allsamples.AddSample("/scratch/fronga/MC_v1.29/ZinvisibleJets_7TeV-madgraph.root","Z nunu",totEventsZnunu,ZnunuCrossSection,false,false,1,kGreen);
489     allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-100To250_7TeV-madgraph.root","QCD",totEventsQCD100to250,QCD100to250CrossSection,false,false,6,kRed);
490     allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-250To500_7TeV-madgraph.root","QCD",totEventsQCD250to500,QCD250to500CrossSection,false,false,6,kRed);
491     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-500To1000_7TeV-madgraph.root","QCD",totEventsQCD500to1000,QCD500to1000CrossSection,false,false,6,kRed);
492     allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-1000_7TeV-madgraph.root","QCD",totEventsQCD1000toInf,QCD1000toInfCrossSection,false,false,6,kRed);
493     allsamples.AddSample("/scratch/fronga/MC_v1.29/LM4_SUSY_sftsht_7TeV-pythia6.root","20xLM4",218536,1.879*20,false,true,7,kMagenta);
494     */
495    
496     //how this works: AddSample("FilePath,"Sample name",Nevents,XS,is_data,is_signal,groupindex
497     // allsamples.AddSample("/scratch/fronga/JSON_160404-163869_v1.29/AllCertified191_v1.29.root","Data",1,1,true,false,0,kBlack);
498     // allsamples.AddSample("/scratch/buchmann/May2011ReRecoAllMayJSON.root","Data",1,1,true,false,0,kBlack);
499    
500     /*
501     allsamples.AddSample("/scratch/buchmann/DataMay2011ReReco_CertifiedReRecoJSON.root","Data",1,1,true,false,0,kBlack);
502     allsamples.AddSample("/scratch/fronga/MC_v1.29/TTJets_TuneZ2_7TeV-madgraph-tauola.root","t#bar{t}",totEventsTTbar,TTbarCrossSection,false,false,1,ttbar_color);
503     allsamples.AddSample("/scratch/fronga/MC_v1.29/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets",totEventsWJets,WJetsCrossSection,false,false,3,wjets_color);
504     allsamples.AddSample("/scratch/fronga/MC_v1.29/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root","Dibosons",totEventsVVJets,VVJetsCrossSection,false,false,4,diboson_color);
505     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_s-channel_7TeV-madgraph_2.root","SingleTop",totEventsSingleTopS,SingleTopSCrossSection,false,false,2,singletop_color);
506     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopT,SingleTopTCrossSection,false,false,2,singletop_color);
507     allsamples.AddSample("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopU,SingleTopUCrossSection,false,false,2,singletop_color);
508     allsamples.AddSample("/scratch/fronga/MC_v1.29/ZinvisibleJets_7TeV-madgraph.root","Z nunu",totEventsZnunu,ZnunuCrossSection,false,false,4,diboson_color);
509     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-100To250_7TeV-madgraph.root","QCD",totEventsQCD100to250,QCD100to250CrossSection,false,false,5,qcd_color);
510     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-250To500_7TeV-madgraph.root","QCD",totEventsQCD250to500,QCD250to500CrossSection,false,false,5,qcd_color);
511     allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-500To1000_7TeV-madgraph.root","QCD",totEventsQCD500to1000,QCD500to1000CrossSection,false,false,5,qcd_color);
512     allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-1000_7TeV-madgraph.root","QCD",totEventsQCD1000toInf,QCD1000toInfCrossSection,false,false,5,qcd_color);
513     allsamples.AddSample("/scratch/fronga/MC_v1.29/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root","ZJets",totEventsZjets,ZJetsCrossSection,false,false,6,dy_color);
514    
515    
516     allsamples.AddSample("/scratch/fronga/MC_v1.29/LM4_SUSY_sftsht_7TeV-pythia6.root","20xLM4",218536,1.879*20,false,true,7,lm_color);
517    
518    
519    
520    
521    
522    
523    
524     // allsamples.AddSample("/scratch/buchmann/DataMay2011ReReco_CertifiedReRecoJSON.root","Data",1,1,true,false,0,kBlack);
525     // allsamples.AddSample("/scratch/buchmann/AllDataJun3_again.root","Data",1,1,true,false,0,kBlack);
526     allsamples.AddSample("/scratch/buchmann/AllData_Jun10___486pb_merged.root","Data",1,1,true,false,0,kBlack);
527     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/TTJets_TuneZ2_7TeV-madgraph-tauola.root","t#bar{t}",totEventsTTbar,TTbarCrossSection,false,false,1,ttbar_color);
528     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root","WJets",totEventsWJets,WJetsCrossSection,false,false,3,wjets_color);
529     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root","Dibosons",totEventsVVJets,VVJetsCrossSection,false,false,4,diboson_color);
530     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/TToBLNu_TuneZ2_s-channel_7TeV-madgraph_2.root","SingleTop",totEventsSingleTopS,SingleTopSCrossSection,false,false,2,singletop_color);
531     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopT,SingleTopTCrossSection,false,false,2,singletop_color);
532     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root","SingleTop",totEventsSingleTopU,SingleTopUCrossSection,false,false,2,singletop_color);
533     // allsamples.AddSample("/scratch/fronga/MC_v1.29/ZinvisibleJets_7TeV-madgraph.root","Z nunu",totEventsZnunu,ZnunuCrossSection,false,false,4,diboson_color);
534     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-100To250_7TeV-madgraph.root","QCD",totEventsQCD100to250,QCD100to250CrossSection,false,false,5,qcd_color);
535     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-250To500_7TeV-madgraph.root","QCD",totEventsQCD250to500,QCD250to500CrossSection,false,false,5,qcd_color);
536     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-500To1000_7TeV-madgraph.root","QCD",totEventsQCD500to1000,QCD500to1000CrossSection,false,false,5,qcd_color);
537     // allsamples.AddSample("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-1000_7TeV-madgraph.root","QCD",totEventsQCD1000toInf,QCD1000toInfCrossSection,false,false,5,qcd_color);
538     allsamples.AddSample("/scratch/buchmann/MC_Spring11_PU_PF/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root","ZJets",totEventsZjets,ZJetsCrossSection,false,false,6,dy_color);
539    
540    
541     // allsamples.AddSample("/scratch/fronga/MC_v1.29/LM4_SUSY_sftsht_7TeV-pythia6.root","20xLM4",218536,1.879*20,false,true,7,lm_color);
542     allsamples.AddSample("/scratch/buchmann/LM4_SUSY_sftsht_7TeV-pythia.root","20xLM4",218536,1.879*20,false,true,7,lm_color);
543    
544    
545    
546    
547    
548    
549    
550    
551    
552    
553    
554    
555    
556    
557    
558     */
559 buchmann 1.2
560    
561     void DoNewExperimentalMetPrediction() {
562     // Do it without Pt correction
563     float MCPeakNoPtCorr=0,MCPeakErrorNoPtCorr=0,DataPeakNoPtCorr=0,DataPeakErrorNoPtCorr=0,MCSigma=0,DataSigma=0;
564     stringstream resultsNoPtCorr;
565     stringstream NoPtCorrdatajzb;
566     stringstream NoPtCorrmcjzb;
567    
568     string originalJZBvariableDATA = PlottingSetup::jzbvariabledata;
569     string originalJZBvariableMC = PlottingSetup::jzbvariablemc;
570    
571     PlottingSetup::jzbvariabledata="jzb[1]";
572     PlottingSetup::jzbvariablemc="jzb[1]";
573    
574     stringstream mcjzbNoPtCorr;
575     find_peaks(MCPeakNoPtCorr,MCPeakErrorNoPtCorr, DataPeakNoPtCorr, DataPeakErrorNoPtCorr,resultsNoPtCorr,true,NoPtCorrdatajzb,NoPtCorrmcjzb);
576     string mcjzb = NoPtCorrmcjzb.str();
577    
578     PlottingSetup::jzbvariabledata=originalJZBvariableDATA;
579     PlottingSetup::jzbvariablemc=originalJZBvariableMC;
580    
581     cout << "Have found a peak in mc at " << MCPeakNoPtCorr << " +/- " << MCPeakErrorNoPtCorr << " ; will use this result." << endl;
582    
583     string CorrectedMet="met[4]+"+any2string(abs(2*(MCPeakNoPtCorr)));
584     if(2*(MCPeakNoPtCorr)<0) CorrectedMet="met[4]-"+any2string(2*(MCPeakNoPtCorr));
585     string CorrectedMetUp="met[4]+"+any2string(abs(2*(MCPeakNoPtCorr+MCPeakErrorNoPtCorr)));
586     if(2*(MCPeakNoPtCorr+MCPeakErrorNoPtCorr)<0) "met[4]-"+any2string(abs(2*(MCPeakNoPtCorr+MCPeakErrorNoPtCorr)));
587     string CorrectedMetDown="met[4]+"+any2string(abs(2*(MCPeakNoPtCorr-MCPeakErrorNoPtCorr)));
588     if(2*(MCPeakNoPtCorr-MCPeakErrorNoPtCorr)<0) "met[4]-"+any2string(abs(2*(MCPeakNoPtCorr+MCPeakErrorNoPtCorr)));
589    
590     do_experimental_zjets_prediction("met[4]",CorrectedMet,"jzb[1]",MCPeakNoPtCorr,"CovertMapping");
591     do_experimental_zjets_prediction("met[4]",CorrectedMetDown,"jzb[1]",(MCPeakNoPtCorr-MCPeakErrorNoPtCorr),"CovertMapping_SysDown");
592     do_experimental_zjets_prediction("met[4]",CorrectedMetUp,"jzb[1]",(MCPeakNoPtCorr+MCPeakErrorNoPtCorr),"CovertMapping_SysUp");
593    
594    
595     }
596    
597     void DoNewExperimentalMetPrediction_WithPtCorrection() {
598     // Do it with Pt correction
599     float MCPeakWithPtCorr=0,MCPeakErrorWithPtCorr=0,DataPeakWithPtCorr=0,DataPeakErrorWithPtCorr=0,MCSigma=0,DataSigma=0;
600     stringstream resultsWithPtCorr;
601     stringstream WithPtCorrdatajzb;
602     stringstream WithPtCorrmcjzb;
603    
604     string originalJZBvariableDATA = PlottingSetup::jzbvariabledata;
605     string originalJZBvariableMC = PlottingSetup::jzbvariablemc;
606    
607     stringstream mcjzbWithPtCorr;
608     // find_peaks(MCPeakWithPtCorr,MCPeakErrorWithPtCorr, DataPeakWithPtCorr, DataPeakErrorWithPtCorr,resultsWithPtCorr,true,WithPtCorrdatajzb,WithPtCorrmcjzb);
609     write_warning(__FUNCTION__,"Hardcoded peak position");MCPeakWithPtCorr=3.58273,MCPeakErrorWithPtCorr=0.208997;
610    
611     string mcjzb = WithPtCorrmcjzb.str();
612    
613     cout << "Have found a peak in mc at " << MCPeakWithPtCorr << " +/- " << MCPeakErrorWithPtCorr << " ; will use this result." << endl;
614    
615     string CorrectedMet="met[4]+"+any2string(abs(2*(MCPeakWithPtCorr)));
616     if(2*(MCPeakWithPtCorr)<0) CorrectedMet="met[4]-"+any2string(2*(MCPeakWithPtCorr));
617     string CorrectedMetUp="met[4]+"+any2string(abs(2*(MCPeakWithPtCorr+MCPeakErrorWithPtCorr)));
618     if(2*(MCPeakWithPtCorr+MCPeakErrorWithPtCorr)<0) "met[4]-"+any2string(abs(2*(MCPeakWithPtCorr+MCPeakErrorWithPtCorr)));
619     string CorrectedMetDown="met[4]+"+any2string(abs(2*(MCPeakWithPtCorr-MCPeakErrorWithPtCorr)));
620     if(2*(MCPeakWithPtCorr-MCPeakErrorWithPtCorr)<0) "met[4]-"+any2string(abs(2*(MCPeakWithPtCorr+MCPeakErrorWithPtCorr)));
621    
622     do_experimental_zjets_prediction("met[4]",CorrectedMet,originalJZBvariableMC,MCPeakWithPtCorr,"CovertMapping_WithPtCorr_AdaptedPredMet");
623     // do_experimental_zjets_prediction("met[4]",CorrectedMetDown,originalJZBvariableMC,(MCPeakWithPtCorr-MCPeakErrorWithPtCorr),"CovertMapping_WithPtCorr_AdaptedPredMet_SysDown");
624     // do_experimental_zjets_prediction("met[4]",CorrectedMetUp,originalJZBvariableMC,(MCPeakWithPtCorr+MCPeakErrorWithPtCorr),"CovertMapping_WithPtCorr_AdaptedPredMet_SysUp");
625    
626     do_experimental_zjets_prediction("met[4]",CorrectedMet+"+0.036*pt",originalJZBvariableMC,MCPeakWithPtCorr,"CovertMapping_WithPtCorr_AdaptedPredMetWithPtDep");
627    
628     // do_experimental_zjets_prediction("met[4]","met[4]",originalJZBvariableMC,MCPeakWithPtCorr,"CovertMapping_WithPtCorr");
629    
630    
631     }
632    
633    
634     void OUTDATED_ExperimentalMetPrediction() {
635    
636     // compute_r_in_out();
637    
638     // DoNewExperimentalMetPrediction();
639     DoNewExperimentalMetPrediction_WithPtCorrection();
640    
641     cout << "About to switch off the weight : " << (const char*) cutWeight << endl;
642     TCut weightbackup=cutWeight;
643     cutWeight="1.0";
644     cout << "Have set cut weight to 1.0 (no PURW)" << endl;
645    
646     // do_experimental_zjets_prediction("met[4]",PlottingSetup::jzbvariablemc,0,"PtCorrectedJZB");
647     /*
648     // Do it without PU re-weighting
649     float MCPeakNoPU=0,MCPeakErrorNoPU=0,DataPeakNoPU=0,DataPeakErrorNoPU=0,MCSigma=0,DataSigma=0;
650     stringstream resultsNoPU;
651     stringstream noPUdatajzb;
652     stringstream noPUmcjzb;
653    
654     stringstream mcjzbnoPU;
655     find_peaks(MCPeakNoPU,MCPeakErrorNoPU, DataPeakNoPU, DataPeakErrorNoPU,resultsNoPU,true,noPUdatajzb,noPUmcjzb);
656     string mcjzb = noPUmcjzb.str();
657     */
658     // do_experimental_zjets_prediction("met[4]","jzb[1]",0,"RAWJZB");
659    
660     // do_experimental_zjets_prediction("met[4]","((jzb[1]))",-0.614166,"PeakCorrectedJZB");/*
661     // do_experimental_zjets_prediction("met[4]",mcjzb,"FullyCorrectedJZB");
662     /*
663     do_experimental_zjets_prediction(tauevents.GetTree(),"met[4]","jzb[1]","tau_RAWJZB");
664     do_experimental_zjets_prediction(tauevents.GetTree(),"met[4]","((jzb[1])- 3.53399 )","tau_PeakCorrectedJZB");
665     do_experimental_zjets_prediction(tauevents.GetTree(),"met[4]","((jzb[1]+0.059979*pt))","tau_PtCorrectedJZB");
666     do_experimental_zjets_prediction(tauevents.GetTree(),"met[4]","((jzb[1]+0.059979*pt)- 3.53399 )","tau_FullyCorrectedJZB");
667     */
668     /**/
669    
670     cutWeight=weightbackup;
671     cout << "Have reset cut weight (PURW reactivated)" << endl;
672     }
673    
674     oid do_experimental_zjets_prediction(string metvariable, string predvariable, string jzbvariable, float epsilonZero, string SaveAsName) {
675     TCanvas *can = new TCanvas("can","can");
676    
677     cout << "Going to use the variable " << metvariable << " to draw met " << endl;
678     cout << "Going to use the variable " << predvariable<< " to draw the met PREDICTION" << endl;
679     cout << "Going to use the variable " << jzbvariable << " to cut on JZB" << endl;
680     int nbins=50;
681     float xmin=0;
682     float xmax=100;
683    
684     string basiccut="pfJetGoodNum40>=3&&abs(mll-91.2)<20&&pt1>20&&pt2>20&&(abs(jzb[1]-"+any2string(epsilonZero)+")>abs("+any2string(epsilonZero)+"))&&((";
685    
686     cout << "The basic cut is " << basiccut << endl;
687     stringstream ObsCommand;
688     ObsCommand << basiccut << jzbvariable << ")>0)&&(id1==id2)&&abs("+metvariable+"-75)>5";
689     stringstream PredCommand;
690     PredCommand << basiccut << jzbvariable << ")<0)&&(id1==id2)&&abs("+predvariable+"-75)>5";
691     stringstream PredCommand3;
692     PredCommand3 << basiccut << jzbvariable << ")>0)&&(id1!=id2)&&abs("+predvariable+"-75)>5";
693     stringstream PredCommand4;
694     PredCommand4 << basiccut << jzbvariable << ")<0)&&(id1!=id2)&&abs("+predvariable+"-75)>5";
695    
696     TH1F *obs = allsamples.Draw("obs",metvariable,nbins,xmin,xmax,"MET [GeV]","events",TCut(ObsCommand.str().c_str()),mc,PlottingSetup::luminosity,allsamples.FindSample("Z_em_DYJetsToLL"));
697     TH1F *pred = allsamples.Draw("pred", predvariable,nbins,xmin,xmax,"MET [GeV]","events",TCut(PredCommand.str().c_str()),mc,PlottingSetup::luminosity,allsamples.FindSample("Z_em_DYJetsToLL"));
698     TH1F *pre3 = allsamples.Draw("pre3",predvariable,nbins,xmin,xmax,"MET [GeV]","events",TCut(PredCommand3.str().c_str()),mc,PlottingSetup::luminosity,allsamples.FindSample("Z_em_DYJetsToLL"));
699     TH1F *pre4 = allsamples.Draw("pre4",predvariable,nbins,xmin,xmax,"MET [GeV]","events",TCut(PredCommand4.str().c_str()),mc,PlottingSetup::luminosity,allsamples.FindSample("Z_em_DYJetsToLL"));
700    
701     obs->SetLineColor(kBlue);
702     obs->SetMarkerSize(0.3);
703     obs->SetMarkerColor(kBlue);
704    
705     pred->Add(pre3);
706     pred->Add(pre4,-1);
707    
708     pred->SetMarkerColor(kRed);
709     pred->SetLineColor(kRed);
710    
711     can->SetLogy(0);
712     obs->Draw("e1");
713     pred->Draw("histo,same");
714     obs->Draw("same,e1");
715    
716     TLegend *leg = new TLegend(0.7,0.7,0.89,0.89);
717     leg->SetFillColor(kWhite);
718     leg->AddEntry(obs,"Observed","lp");
719     leg->AddEntry(pred,"Predicted","lp");
720     leg->Draw("same");
721    
722     stringstream SaveAsName1;
723     SaveAsName1 << "MetPredPlots/" << SaveAsName;
724     CompleteSave(can,SaveAsName1.str());
725    
726     can->SetLogy(0);
727     TH1F *ratio = (TH1F*)obs->Clone("ratio");
728     ratio->SetTitle("obs/pred");
729     ratio->SetName("ratio");
730     ratio->Divide(pred);
731     ratio->GetYaxis()->SetTitle("obs/pred");
732     ratio->GetYaxis()->CenterTitle();
733     ratio->GetYaxis()->SetRangeUser(0,5);
734     TF1 *func = new TF1("func","[0]+[1]*x",10,100);
735     ratio->Fit("func","R");
736     ratio->Draw();
737     TLine *line = new TLine(0.0,1.0,100,1.0);
738     line->SetLineColor(kBlue);
739     line->SetLineWidth(2);
740     line->Draw("same");
741    
742     stringstream SaveAsName2;
743     SaveAsName2 << "MetPredPlots/" << SaveAsName << "_ratio";
744     CompleteSave(can,SaveAsName2.str());
745    
746     delete line;
747     delete ratio;
748     delete pred;
749     delete pre3;
750     delete pre4;
751     delete obs;
752     delete can;
753     }