ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/cbrown/AnalysisFramework/Mostas/Mostas.C
Revision: 1.1
Committed: Wed Jun 22 11:38:19 2011 UTC (13 years, 10 months ago) by buchmann
Content type: text/plain
Branch: MAIN
CVS Tags: cbaf_4_98ifb_paper, cbaf_4p7ifb, Honeypot, cbaf_2p1ifb, HEAD
Log Message:
Initial commit of Mostas script

File Contents

# User Rev Content
1 buchmann 1.1 /*
2    
3     load this in root and execute it with .x Mostas.C+
4    
5     this script is also found as void smeared_simple_2011Jun10DCS_PF()
6    
7     */
8     #include <iostream>
9     #include <string>
10     #include <vector>
11     #include <sstream>
12     #include <iomanip>
13     #include <algorithm>
14    
15     #include "TFile.h"
16     #include "TLine.h"
17     #include "TCanvas.h"
18     #include "TCut.h"
19     #include "TTree.h"
20     #include "TH1F.h"
21     #include "TH2F.h"
22     #include "TStyle.h"
23     #include "TROOT.h"
24     #include "TF1.h"
25     #include "TMath.h"
26     #include "TGaxis.h"
27     #include "TLegend.h"
28     #include "THStack.h"
29     #include "TPaveText.h"
30     #include "TColor.h"
31     #include "TGraphAsymmErrors.h"
32     #include "TRandom.h"
33    
34     #include "setTDRStyle.C"
35    
36     using namespace std;
37    
38     #define nSamples 20
39    
40     TFile *file[nSamples];
41     TTree *events[nSamples];
42     float eventsWeight[nSamples]; // for MC not yet automatic available, for data eventWeight = 1
43     vector<TFile*> files;
44     vector<string> filenames;
45    
46     vector<string> histTitleList;
47    
48     TH1F * sumJetPt[nSamples];
49    
50     TH1F * ZPt[nSamples];
51    
52     float lumi=648.0;//486.0;//
53    
54     /*
55    
56     2010 : 34.1 /pb, 23087 Z events with abs(mll-91.2)<20&&pfJetGoodNum>=0
57     2011: 22894 Z events with abs(mll-91.2)<20&&pfJetGoodNum>=0 so --> 33.8 /pb :-)
58    
59    
60     */
61    
62    
63    
64     TH1F * J1Pt_ossf[nSamples];
65     TH1F * J1Eta_ossf[nSamples];
66     TH1F * ZPt_ossf_incl[nSamples];
67     TH1F * leptonPt[nSamples];
68     TH1F * leptonEta[nSamples];
69     TH1F * mll[nSamples];
70     TH1F * mll_ossf[nSamples];
71     TH1F * mll_osof[nSamples];
72    
73     TH1F * jzb_ossf[nSamples];
74     TH1F * jzb_ossfee[nSamples];
75     TH1F * jzb_ossfmm[nSamples];
76     TH1F * JPt_ossf[nSamples];
77     TH1F * ZPt_ossf[nSamples];
78     TH1F * nJets[nSamples];
79     TH1F * nJets285[nSamples];
80     TH1F * nJets315[nSamples];
81     TH1F * jzb_osof[nSamples];
82     TH1F * mll_ee[nSamples];
83     TH1F * mll_mumu[nSamples];
84    
85    
86     TH1F * jzbprime_ossf[nSamples];
87     TH1F * jzbprime_osof[nSamples];
88    
89     TH1F * sjzbprime_ossf[nSamples];
90    
91     TH1F * jzbprime_ossfee[nSamples];
92     TH1F * jzbprime_ossfmm[nSamples];
93    
94     TH1F * jzbprime_ossf_JZBP[nSamples];
95     TH1F * jzbprime_ossf_JZBN[nSamples];
96     TH1F * jzbprime_osof_JZBP[nSamples];
97     TH1F * jzbprime_osof_JZBN[nSamples];
98    
99     TH1F * jzbprime_ossf_JZBPfin[nSamples];
100     TH1F * jzbprime_ossf_JZBNfin[nSamples];
101     TH1F * jzbprime_osof_JZBPfin[nSamples];
102     TH1F * jzbprime_osof_JZBNfin[nSamples];
103    
104     TH1F * jzbprime_ossf_JZBP_MisCalibP[nSamples];
105     TH1F * jzbprime_ossf_JZBN_MisCalibP[nSamples];
106     TH1F * jzbprime_osof_JZBP_MisCalibP[nSamples];
107     TH1F * jzbprime_osof_JZBN_MisCalibP[nSamples];
108    
109     TH1F * jzbprime_ossf_JZBP_MisCalibN[nSamples];
110     TH1F * jzbprime_ossf_JZBN_MisCalibN[nSamples];
111     TH1F * jzbprime_osof_JZBP_MisCalibN[nSamples];
112     TH1F * jzbprime_osof_JZBN_MisCalibN[nSamples];
113    
114     // not auto created
115     TH1F * jzb_subtracted[nSamples];
116     TH1F * jzb_subtractedee[nSamples];
117     TH1F * jzb_subtractedmm[nSamples];
118    
119     TH1F * Bpred[nSamples];
120     TH1F * Bpredfin[nSamples];
121     TH1F * BpredSysP[nSamples];
122     TH1F * BpredSysN[nSamples];
123     //TH1F * JZBratio[nSamples];
124     TGraphAsymmErrors *JZBratiofin[nSamples];
125     TGraphAsymmErrors * JZBratio[nSamples];
126     TH1F * JZBpull[nSamples];
127    
128     TH1F * jzb_ossf_rebinned[nSamples];
129     TH1F * jzb_osof_rebinned[nSamples];
130    
131    
132    
133     TF1 *f1MC;
134     TF1 *f1Data;
135     TF1 *pol0[nSamples];
136     TF1 *cb[nSamples];
137     TF1 *cbP[nSamples];
138     TF1 *cbN[nSamples];
139    
140     int canvasCounter;
141    
142    
143     Int_t nice_blue ;
144     Int_t nice_green ;
145     Int_t nice_red ;
146     Int_t nice_pink ;
147     Int_t nice_orange;//#FA9624
148     Int_t nice_black ;//black with a bit of purple :-)
149    
150    
151    
152    
153     TLegend *globalDataMCLeg;
154     TLegend *globalMCLeg;
155    
156     TPaveText *paveTitle;
157     TPaveText *paveTitleMCZjets;
158     TPaveText *paveTitleMC;
159     TPaveText *ij2jzbpPave;
160     TPaveText *ij2jzbnPave;
161    
162    
163     float statErrorN(float x){return x - 0.5*TMath::ChisquareQuantile(0.3173/2,2*x);}
164     float statErrorP(float x){return 0.5*TMath::ChisquareQuantile(1-0.3173/2,2*(x+1))-x;}
165     float lowLimit(float a, float x){return 0.5*TMath::ChisquareQuantile(a,2*x);}
166     float highLimit(float a,float x){return 0.5*TMath::ChisquareQuantile(1-a,2*(x+1));}
167    
168     THStack *mc_stack(TH1F** histos);
169     THStack *mc_stack(TH1F** histos,int rebin);
170     THStack *mc_normalized_stack(TH1F** histos);
171     THStack *mc_normalizedtodata_stack(TH1F** histos);
172     THStack *mc_normalizedtodata_stack(TH1F** histos, int rebin); //blublu
173     TH1F *mc_th1f(TH1F** histos);
174     TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id);
175     TGraphAsymmErrors *histRatiofin(TH1F *h1,TH1F *h2, int id);
176    
177     string fileNamePreffix;
178     void setYMax(TH1F*hist,float x);
179     void setYMin(TH1F*hist,float x);
180    
181     // sample - leptons - mllcut - nJets - variable
182     enum sampleIndex {data, allB, zjets, ttbar, lm4, zjets_un, wjets, znunu, vvjets, singletops, singletopt, singletopu, qcd100to250, qcd250to500, qcd500to1000, qcd1000toInf, BnoZ,BwithZ,BandS};
183    
184     vector<int> allBList;
185     vector<TH1F**> allTH1F;
186    
187     template<class T> std::string any2string(T i);
188    
189     TCanvas *can[100];
190     TLegend *leg[100];
191     TPaveText *paveText[100];
192     THStack *hs[100];
193     THStack *hsErrorBand[100];
194     TLine *line[100][100];
195     TF1 *f1line[100][100];
196    
197     void saveCanvas(int);
198     void saveCanvas(int, string filename);
199     float histIntegral(TH1F *hist, float minMet);
200     float histIntegralAndError(TH1F *hist, float minMet, float & error);
201     float histIntegral(TH1F *hist, float minMet, float maxMet);
202     float histIntegralAndError(TH1F *hist, float minMet, float maxMet, float & error);
203     float computeRatioError(float a, float da, float b, float db);
204     float computeProductError(float a, float da, float b, float db);
205    
206     Double_t InvCrystalBall(Double_t *x,Double_t *par);
207     Double_t InvCrystalBallP(Double_t *x,Double_t *par);
208     Double_t InvCrystalBallN(Double_t *x,Double_t *par);
209     Double_t CrystalBall(Double_t *x,Double_t *par);
210    
211    
212     TH1F *agreementHist(TH1F *ratioPlot);
213     TH1F *agreementHist(TH1F *obs, TH1F *pred);
214    
215     //void compute_Gaussian(TH1F *hist, TF1 *g, float &err);
216     //float get_peak(TH1F *histo, float &err);
217     float get_peak(TH1F *histo, float &err, float &sigma, TF1*fitFunc, float lowlimit, float highlimit);
218    
219     Double_t GausRandom(Double_t mu, Double_t sigma) {
220     return gRandom->Gaus(mu,sigma); //real deal
221     //return mu;//debugging : no smearing
222     }
223    
224     class sample{
225     public:
226     sample(){;}
227     sample(string fileName,float w, int index){fileName_ = fileName;w_=w; fp_ = new TFile(fileName_.c_str(),"OPEN"); tree_ = (TTree*)fp_->Get("PFevents");weight_=any2string(w_).c_str();index_=index;}
228     string fileName_;
229     TFile *fp_;
230     TTree *tree_;
231     float w_;
232     int index_;
233     TCut weight_;
234     };
235    
236     //const string JZBtype = "sumJetPt[1]-pt";
237     const string JZBtype = "jzb[1]";
238     //const string JZBtype = "jzb[2]";
239     const string METtype = "met[4]";
240     // nEvent.jzb[2] = recoil.Pt() - (s1+s2).Pt(); // to be used recoil met (recoilpt[0])
241     // nEvent.jzb[1] = pfNoCutsJetVector.Pt() - (s1+s2).Pt(); // to be used with pfMET
242     // nEvent.jzb[4] = tcNoCutsJetVector.Pt() - (s1+s2).Pt(); // to be used with tcMET
243    
244    
245     float fin_binning[]={0,5,10,20,35,50,100,350};
246    
247     void createHisto(int sampleIndex)
248     {
249     string histTitle;
250    
251    
252     int jetPtBins = 40;
253     float jetMin = 0;
254     float jetMax = 200;
255    
256     int JZBbins=1400;
257     // int JZBbins=350;
258     int JZBcorserbins=280; // nJets>=2
259     // int JZBcorserbins=350; // nJets>=2
260     int absJZBcorserbins=JZBcorserbins/2; // nJets>=2
261     int JZBmin=-700;
262     int JZBmax=700;
263    
264     //float absJZBbins[] = {0,5,10,20,35,55,80,110,140};
265     //int absJZBbinsNumber = sizeof(absJZBbins)/4-1;
266    
267     int mllBins =100;
268     float mllMin = 20.;
269     float mllMax =220.;
270    
271    
272     nice_blue = TColor::GetColor("#2E9AFE");
273     nice_green = TColor::GetColor("#81f781");
274     nice_red = TColor::GetColor("#F78181");
275     nice_pink = TColor::GetColor("#F781BE");
276     nice_orange= TColor::GetColor("#F7BE81");//#FA9624
277     nice_black = TColor::GetColor("#2A0A1B");//black with a bit of purple :-)
278    
279    
280    
281     histTitle = "sumJetPt_";
282     histTitleList.push_back(histTitle);
283     histTitle = histTitle+any2string(sampleIndex);
284     sumJetPt[sampleIndex] = new TH1F(histTitle.c_str(),"; vector sum jet p_{T} (GeV) ; events",jetPtBins,jetMin,jetMax);
285     sumJetPt[sampleIndex]->Sumw2();
286     if(sampleIndex==allB)allTH1F.push_back(sumJetPt);
287    
288    
289     histTitle = "ZPt_";
290     histTitleList.push_back(histTitle);
291     histTitle = histTitle+any2string(sampleIndex);
292     ZPt[sampleIndex] = new TH1F(histTitle.c_str(),"; Z p_{T} (GeV) ; events",jetPtBins,jetMin,jetMax);
293     ZPt[sampleIndex]->Sumw2();
294     if(sampleIndex==allB)allTH1F.push_back(ZPt);
295    
296     histTitle = "J1Pt_ossf_";
297     histTitleList.push_back(histTitle);
298     histTitle = histTitle+any2string(sampleIndex);
299     J1Pt_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; leading jet Pt (GeV) ; events",jetPtBins,jetMin,jetMax);
300     J1Pt_ossf[sampleIndex]->Sumw2();
301     if(sampleIndex==allB)allTH1F.push_back(J1Pt_ossf);
302    
303    
304     histTitle = "J1Eta_ossf_";
305     histTitleList.push_back(histTitle);
306     histTitle = histTitle+any2string(sampleIndex);
307     J1Eta_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; leading jet #eta ; events",40,-5,5);
308     J1Eta_ossf[sampleIndex]->Sumw2();
309     if(sampleIndex==allB)allTH1F.push_back(J1Eta_ossf);
310    
311     histTitle = "ZPt_ossf_incl_";
312     histTitleList.push_back(histTitle);
313     histTitle = histTitle+any2string(sampleIndex);
314     ZPt_ossf_incl[sampleIndex] = new TH1F(histTitle.c_str(),"; Z p_{T} (GeV) ; events",jetPtBins,jetMin,jetMax);
315     ZPt_ossf_incl[sampleIndex]->Sumw2();
316     if(sampleIndex==allB)allTH1F.push_back(ZPt_ossf_incl);
317    
318    
319     histTitle = "leptonPt_";
320     histTitleList.push_back(histTitle);
321     histTitle = histTitle+any2string(sampleIndex);
322     leptonPt[sampleIndex] = new TH1F(histTitle.c_str(),"; p_{T} (GeV) ; events",100,0,200);
323     leptonPt[sampleIndex]->Sumw2();
324     if(sampleIndex==allB)allTH1F.push_back(leptonPt);
325    
326     histTitle = "leptonEta_";
327     histTitleList.push_back(histTitle);
328     histTitle = histTitle+any2string(sampleIndex);
329     leptonEta[sampleIndex] = new TH1F(histTitle.c_str(),"; #eta ; events",60,-3.0,3.0);
330     leptonEta[sampleIndex]->Sumw2();
331     if(sampleIndex==allB)allTH1F.push_back(leptonEta);
332    
333     histTitle = "mll_";
334     histTitleList.push_back(histTitle);
335     histTitle = histTitle+any2string(sampleIndex);
336     mll[sampleIndex] = new TH1F(histTitle.c_str(),"; invariant mass (GeV^{2}) ; events",200,50,250);
337     mll[sampleIndex]->Sumw2();
338     if(sampleIndex==allB)allTH1F.push_back(mll);
339    
340    
341    
342     histTitle = "mll_osof_";
343     histTitleList.push_back(histTitle);
344     histTitle = histTitle+any2string(sampleIndex);
345     mll_osof[sampleIndex] = new TH1F(histTitle.c_str(),"; inv. mass (GeV) ; events",mllBins,mllMin,mllMax);
346     mll_osof[sampleIndex]->Sumw2();
347     if(sampleIndex==allB)allTH1F.push_back(mll_osof);
348    
349     histTitle = "mll_ossf_";
350     histTitleList.push_back(histTitle);
351     histTitle = histTitle+any2string(sampleIndex);
352     mll_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; inv. mass (GeV) ; events",mllBins,mllMin,mllMax);
353     mll_ossf[sampleIndex]->Sumw2();
354     if(sampleIndex==allB)allTH1F.push_back(mll_ossf);
355    
356    
357     histTitle = "JPt_ossf_";
358     histTitleList.push_back(histTitle);
359     histTitle = histTitle+any2string(sampleIndex);
360     JPt_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; leading jet Pt (GeV) ; events",jetPtBins,jetMin,jetMax);
361     JPt_ossf[sampleIndex]->Sumw2();
362     if(sampleIndex==allB)allTH1F.push_back(JPt_ossf);
363    
364     histTitle = "ZPt_ossf_";
365     histTitleList.push_back(histTitle);
366     histTitle = histTitle+any2string(sampleIndex);
367     ZPt_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; Z p_{T} (GeV) ; events",jetPtBins,jetMin,jetMax);
368     ZPt_ossf[sampleIndex]->Sumw2();
369     if(sampleIndex==allB)allTH1F.push_back(ZPt_ossf);
370    
371     histTitle = "nJets_";
372     histTitleList.push_back(histTitle);
373     histTitle = histTitle+any2string(sampleIndex);
374     nJets[sampleIndex] = new TH1F(histTitle.c_str(),"; nJets ; events",15,0,15);
375     nJets[sampleIndex]->Sumw2();
376     if(sampleIndex==allB)allTH1F.push_back(nJets);
377    
378     histTitle = "nJets285_";
379     histTitleList.push_back(histTitle);
380     histTitle = histTitle+any2string(sampleIndex);
381     nJets285[sampleIndex] = new TH1F(histTitle.c_str(),"; nJets (>28.5 GeV) ; events",15,0,15);
382     nJets285[sampleIndex]->Sumw2();
383     if(sampleIndex==allB)allTH1F.push_back(nJets285);
384    
385     histTitle = "nJets315_";
386     histTitleList.push_back(histTitle);
387     histTitle = histTitle+any2string(sampleIndex);
388     nJets315[sampleIndex] = new TH1F(histTitle.c_str(),"; nJets (>31.5 GeV) ; events",15,0,15);
389     nJets315[sampleIndex]->Sumw2();
390     if(sampleIndex==allB)allTH1F.push_back(nJets315);
391    
392     histTitle = "jzb_ossf_";
393     histTitleList.push_back(histTitle);
394     histTitle = histTitle+any2string(sampleIndex);
395     jzb_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBbins,JZBmin,JZBmax);
396     jzb_ossf[sampleIndex]->Sumw2();
397     if(sampleIndex==allB)allTH1F.push_back(jzb_ossf);
398    
399     histTitle = "jzb_ossf_ee";
400     histTitleList.push_back(histTitle);
401     histTitle = histTitle+any2string(sampleIndex);
402     jzb_ossfee[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBbins,JZBmin,JZBmax);
403     jzb_ossfee[sampleIndex]->Sumw2();
404     if(sampleIndex==allB)allTH1F.push_back(jzb_ossfee);
405    
406     histTitle = "jzb_ossf_mm";
407     histTitleList.push_back(histTitle);
408     histTitle = histTitle+any2string(sampleIndex);
409     jzb_ossfmm[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBbins,JZBmin,JZBmax);
410     jzb_ossfmm[sampleIndex]->Sumw2();
411     if(sampleIndex==allB)allTH1F.push_back(jzb_ossfmm);
412    
413     /* histTitle = "jzb_ossffin";
414     histTitleList.push_back(histTitle);
415     histTitle = histTitle+any2string(sampleIndex);
416     jzb_ossffin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",8,{0,5,10,20,35,50,75,100,350});
417     jzb_ossffin[sampleIndex]->Sumw2();
418     if(sampleIndex==allB)allTH1F.push_back(jzb_ossffin);
419    
420     histTitle = "jzb_osoffin";
421     histTitleList.push_back(histTitle);
422     histTitle = histTitle+any2string(sampleIndex);
423     jzb_ossffin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",8,{0,5,10,20,35,50,75,100,350});
424     jzb_ossffin[sampleIndex]->Sumw2();
425     if(sampleIndex==allB)allTH1F.push_back(jzb_ossffin);
426     */
427     histTitle = "jzb_osof_";
428     histTitleList.push_back(histTitle);
429     histTitle = histTitle+any2string(sampleIndex);
430     jzb_osof[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBbins,JZBmin,JZBmax);
431     jzb_osof[sampleIndex]->Sumw2();
432     if(sampleIndex==allB)allTH1F.push_back(jzb_osof);
433    
434     histTitle = "jzbprime_ossf_";
435     histTitleList.push_back(histTitle);
436     histTitle = histTitle+any2string(sampleIndex);
437     jzbprime_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBcorserbins,JZBmin,JZBmax);
438     jzbprime_ossf[sampleIndex]->Sumw2();
439     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf);
440    
441     histTitle = "sjzbprime_ossf_";
442     histTitleList.push_back(histTitle);
443     histTitle = histTitle+any2string(sampleIndex);
444     sjzbprime_ossf[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBcorserbins,JZBmin,JZBmax);
445     sjzbprime_ossf[sampleIndex]->Sumw2();
446     if(sampleIndex==allB)allTH1F.push_back(sjzbprime_ossf);
447    
448    
449     histTitle = "jzbprime_osof_";
450     histTitleList.push_back(histTitle);
451     histTitle = histTitle+any2string(sampleIndex);
452     jzbprime_osof[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",JZBcorserbins,JZBmin,JZBmax);
453     jzbprime_osof[sampleIndex]->Sumw2();
454     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof);
455    
456    
457     histTitle = "jzbprime_ossf_JZBP_";
458     histTitleList.push_back(histTitle);
459     histTitle = histTitle+any2string(sampleIndex);
460     jzbprime_ossf_JZBP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
461     jzbprime_ossf_JZBP[sampleIndex]->Sumw2();
462     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBP);
463    
464     histTitle = "jzbprime_ossf_JZBN_";
465     histTitleList.push_back(histTitle);
466     histTitle = histTitle+any2string(sampleIndex);
467     jzbprime_ossf_JZBN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
468     jzbprime_ossf_JZBN[sampleIndex]->Sumw2();
469     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBN);
470    
471     histTitle = "jzbprime_osof_JZBP_";
472     histTitleList.push_back(histTitle);
473     histTitle = histTitle+any2string(sampleIndex);
474     jzbprime_osof_JZBP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
475     jzbprime_osof_JZBP[sampleIndex]->Sumw2();
476     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBP);
477    
478     histTitle = "jzbprime_osof_JZBN_";
479     histTitleList.push_back(histTitle);
480     histTitle = histTitle+any2string(sampleIndex);
481     jzbprime_osof_JZBN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
482     jzbprime_osof_JZBN[sampleIndex]->Sumw2();
483     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBN);
484    
485     histTitle = "jzbprime_ossf_JZBPfin_";
486     histTitleList.push_back(histTitle);
487     histTitle = histTitle+any2string(sampleIndex);
488     jzbprime_ossf_JZBPfin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",7,fin_binning);
489     jzbprime_ossf_JZBPfin[sampleIndex]->Sumw2();
490     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBPfin);
491    
492     histTitle = "jzbprime_ossf_JZBNfin_";
493     histTitleList.push_back(histTitle);
494     histTitle = histTitle+any2string(sampleIndex);
495     jzbprime_ossf_JZBNfin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",7,fin_binning);
496     jzbprime_ossf_JZBNfin[sampleIndex]->Sumw2();
497     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBNfin);
498    
499     histTitle = "jzbprime_osof_JZBPfin_";
500     histTitleList.push_back(histTitle);
501     histTitle = histTitle+any2string(sampleIndex);
502     jzbprime_osof_JZBPfin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",7,fin_binning);
503     jzbprime_osof_JZBPfin[sampleIndex]->Sumw2();
504     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBPfin);
505    
506     histTitle = "jzbprime_osof_JZBNfin_";
507     histTitleList.push_back(histTitle);
508     histTitle = histTitle+any2string(sampleIndex);
509     jzbprime_osof_JZBNfin[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",7,fin_binning);
510     jzbprime_osof_JZBNfin[sampleIndex]->Sumw2();
511     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBNfin);
512    
513    
514     histTitle = "jzbprime_ossf_JZBP_MisCalibP_";
515     histTitleList.push_back(histTitle);
516     histTitle = histTitle+any2string(sampleIndex);
517     jzbprime_ossf_JZBP_MisCalibP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
518     jzbprime_ossf_JZBP_MisCalibP[sampleIndex]->Sumw2();
519     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBP_MisCalibP);
520    
521     histTitle = "jzbprime_ossf_JZBN_MisCalibP_";
522     histTitleList.push_back(histTitle);
523     histTitle = histTitle+any2string(sampleIndex);
524     jzbprime_ossf_JZBN_MisCalibP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
525     jzbprime_ossf_JZBN_MisCalibP[sampleIndex]->Sumw2();
526     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBN_MisCalibP);
527    
528     histTitle = "jzbprime_osof_JZBP_MisCalibP_";
529     histTitleList.push_back(histTitle);
530     histTitle = histTitle+any2string(sampleIndex);
531     jzbprime_osof_JZBP_MisCalibP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
532     jzbprime_osof_JZBP_MisCalibP[sampleIndex]->Sumw2();
533     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBP_MisCalibP);
534    
535     histTitle = "jzbprime_osof_JZBN_MisCalibP_";
536     histTitleList.push_back(histTitle);
537     histTitle = histTitle+any2string(sampleIndex);
538     jzbprime_osof_JZBN_MisCalibP[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
539     jzbprime_osof_JZBN_MisCalibP[sampleIndex]->Sumw2();
540     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBN_MisCalibP);
541    
542    
543     histTitle = "jzbprime_ossf_JZBP_MisCalibN_";
544     histTitleList.push_back(histTitle);
545     histTitle = histTitle+any2string(sampleIndex);
546     jzbprime_ossf_JZBP_MisCalibN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
547     jzbprime_ossf_JZBP_MisCalibN[sampleIndex]->Sumw2();
548     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBP_MisCalibN);
549    
550     histTitle = "jzbprime_ossf_JZBN_MisCalibN_";
551     histTitleList.push_back(histTitle);
552     histTitle = histTitle+any2string(sampleIndex);
553     jzbprime_ossf_JZBN_MisCalibN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
554     jzbprime_ossf_JZBN_MisCalibN[sampleIndex]->Sumw2();
555     if(sampleIndex==allB)allTH1F.push_back(jzbprime_ossf_JZBN_MisCalibN);
556    
557     histTitle = "jzbprime_osof_JZBP_MisCalibN_";
558     histTitleList.push_back(histTitle);
559     histTitle = histTitle+any2string(sampleIndex);
560     jzbprime_osof_JZBP_MisCalibN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
561     jzbprime_osof_JZBP_MisCalibN[sampleIndex]->Sumw2();
562     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBP_MisCalibN);
563    
564     histTitle = "jzbprime_osof_JZBN_MisCalibN_";
565     histTitleList.push_back(histTitle);
566     histTitle = histTitle+any2string(sampleIndex);
567     jzbprime_osof_JZBN_MisCalibN[sampleIndex] = new TH1F(histTitle.c_str(),"; JZB (GeV) ; events",absJZBcorserbins,0,JZBmax);
568     jzbprime_osof_JZBN_MisCalibN[sampleIndex]->Sumw2();
569     if(sampleIndex==allB)allTH1F.push_back(jzbprime_osof_JZBN_MisCalibN);
570    
571    
572     histTitle = "mll_ee_";
573     histTitleList.push_back(histTitle);
574     histTitle = histTitle+any2string(sampleIndex);
575     mll_ee[sampleIndex] = new TH1F(histTitle.c_str(),"; inv. mass (GeV) ; events",mllBins,mllMin,mllMax);
576     mll_ee[sampleIndex]->Sumw2();
577     if(sampleIndex==allB)allTH1F.push_back(mll_ee);
578    
579     histTitle = "mll_mumu_";
580     histTitleList.push_back(histTitle);
581     histTitle = histTitle+any2string(sampleIndex);
582     mll_mumu[sampleIndex] = new TH1F(histTitle.c_str(),"; inv. mass (GeV) ; events",mllBins,mllMin,mllMax);
583     mll_mumu[sampleIndex]->Sumw2();
584     if(sampleIndex==allB)allTH1F.push_back(mll_mumu);
585    
586     }
587     void newCanvas(int&);
588    
589     THStack* errorBand(TH1F* h1, TH1F* h1N, TH1F* h1P ) {
590    
591     THStack *stack = new THStack(TString("stack")+h1->GetName(),h1->GetTitle());
592     h1N->SetLineColor(kWhite);
593     stack->Add(h1N);
594     TH1F* temp = (TH1F*)h1P->Clone();
595     temp->Add(h1N,-1);
596     temp->SetFillColor(kGray);
597     temp->SetLineColor(kWhite);
598     stack->Add(temp);
599     return stack;
600    
601     }
602    
603     void Mostas()
604     //int main()
605     {
606     setTDRStyle();
607     string filename="";
608     string path = "/shome/theofil/JZBTreeData/";
609     gStyle->SetOptTitle(0);
610     gStyle->SetOptFit(0);
611     // gStyle->SetOptStat(1111111);
612     gStyle->SetOptStat(0);
613     TGaxis::SetMaxDigits(3);
614    
615    
616     //float lumi=11.68; // normalize MC to lumi
617    
618     float ZJetsCrossSection = 3048.0; //NNLO----------------------------------
619     float TTbarCrossSection = 165.0;//(NLO) ---- 165.0; // approx. NNLO-----
620     float WJetsCrossSection = 31314.0;//NNLO-------3.131e4; //NNLO------------
621     float ZnunuCrossSection = 5760.0;//NNLO -------4.5e+3; //(LO);------------
622     float SingleTopSCrossSection = 1.49; // NLO;----------------------------------
623     float SingleTopTCrossSection = 20.540; // NLO;--------------------------------
624     float SingleTopUCrossSection = 10.6; // NLO;----------------------------------
625     float VVJetsCrossSection = 4.8; // LO;-------------------------------------
626     float LM4CrossSection = 2.537*20.0; // k*LO
627     /*
628     float QCD100to250CrossSection = 7e+6 ; // LO
629     float QCD250to500CrossSection = 1.71e+5 ; // LO
630     float QCD500to1000CrossSection = 5.2e+3 ; // LO
631     float QCD1000toInfCrossSection = 8.3e+1 ; // LO
632     */
633     //the following numbers are from the MadGraphStandardModel210Summary as linked on the GeneratorProduction2011 page
634     // float QCD50to100CrossSection=30000000; // not used
635     float QCD100to250CrossSection=7000000;
636     float QCD250to500CrossSection=171000;
637     float QCD500to1000CrossSection=5200;
638     float QCD1000toInfCrossSection=83.0;
639    
640     // float totEventsQCD50to100=207418.0; // not used
641     float totEventsQCD100to250=638792.0;
642     float totEventsQCD250to500=344454.0;
643     float totEventsQCD500to1000=10e10;//unknown but now suppressed.
644     float totEventsQCD1000toInf=163185.0;
645    
646     float totEventsZjets = 2313911.0;//2743142.0; // DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola
647     float totEventsTTbar = 1161621.0;//1144028.0; // PabloV13/TTJets_TuneZ2_7TeV-madgraph-tauola
648     float totEventsWJets = 15010237.0;//14350756.0; /// WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root
649     float totEventsZnunu = 2106977.0;//2167964.0; //ZinvisibleJets_7TeV-madgraph.root
650     float totEventsVVJets = 959076.0;//509072.0;
651     float totEventsSingleTopS = 493868.0;//489472.0;
652     float totEventsSingleTopT = 475460.0;//477610.0;
653     float totEventsSingleTopU = 489417.0;//477599.0;
654     /*
655     float totEventsQCD100to250 = 10820757.0;
656     float totEventsQCD250to500 = 4642821.;
657     float totEventsQCD500to1000 = 3152612.;
658     float totEventsQCD1000toInf = 479830.;
659     */
660    
661     cout<<"#--- starting JZB analysis"<<endl;
662    
663    
664     // sample Data("/scratch/buchmann/ntuples2011/data/DataMay201120110509_nightlybuild_2217-AllData.root",1.0,data);//DCS only ...
665     // sample Data("/scratch/buchmann/AllCertified191_v2.root",1.0,data);//191 /pb, duplicates removed, new jet definitions
666     //sample Data("/scratch/fronga/noJSON_v1.29/AllDCS_v1.29.root",1.0,data);//DCS, with Kostas flavor, new code rev1.29
667     // sample Data("/scratch/fronga/JSON_160404-163869_v1.29/AllCertified191_v1.29.root",1.0,data);//DCS, with Kostas flavor, new code rev1.29
668     // sample Data("/scratch/buchmann/Data20110603_GoldenJSON2/AllData_v1p30.root",1.0,data);//336 /pb (Golden JSON, June 3)
669     // sample Data("/scratch/buchmann/AllData_Jun10___486pb_merged.root",1.0,data);//Re-Recoed Data (2011), ~191 /pb or 231 /pb? who knows....
670     // sample Data("/scratch/buchmann/AllData_Jun10___486pb_MoreTriggers3.root",1.0,data);//all data acc2 json jun10
671     // sample Data("/scratch/buchmann/AllData_Jun10___486pb_MoreTriggers3__MAYRERECOONLY.root",1.0,data);//re-reco (may) only!
672    
673     sample Data("/scratch/buchmann/AllData_Jun10___486pb_MoreTriggers4_DCS.root",1.0,data);//DCS, June 10 (648/pb)
674     sample Zjets("/scratch/fronga/MC_v1.29/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root",(lumi/totEventsZjets)*ZJetsCrossSection,zjets);
675     sample Zjets_un("/scratch/fronga/MC_v1.29/DYJetsToLL_TuneZ2_M-50_7TeV-madgraph-tauola.root",1.0,zjets_un); // emulate full MC lumi ~900 pb-1 for Z+jets
676     sample TTbar("/scratch/fronga/MC_v1.29/TTJets_TuneZ2_7TeV-madgraph-tauola.root",((lumi/totEventsTTbar)*TTbarCrossSection),ttbar);
677     sample WJets("/scratch/fronga/MC_v1.29/WJetsToLNu_TuneZ2_7TeV-madgraph-tauola.root",((lumi/totEventsWJets)*WJetsCrossSection),wjets);
678     sample Znunu("/scratch/fronga/MC_v1.29/ZinvisibleJets_7TeV-madgraph.root",((lumi/totEventsZnunu)*ZnunuCrossSection),znunu);
679     sample VVJets("/scratch/fronga/MC_v1.29/VVJetsTo4L_TuneD6T_7TeV-madgraph-tauola.root",((lumi/totEventsVVJets)*VVJetsCrossSection),vvjets);
680     sample SingleTopS("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_s-channel_7TeV-madgraph_2.root",((lumi/totEventsSingleTopS)*SingleTopSCrossSection),singletops);
681     sample SingleTopT("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_t-channel_7TeV-madgraph.root ",((lumi/totEventsSingleTopT)*SingleTopTCrossSection),singletopt);
682     sample SingleTopU("/scratch/fronga/MC_v1.29/TToBLNu_TuneZ2_tW-channel_7TeV-madgraph.root",((lumi/totEventsSingleTopU)*SingleTopUCrossSection),singletopu);
683     sample QCD100to250("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-100To250_7TeV-madgraph.root",((lumi/totEventsQCD100to250)*QCD100to250CrossSection),qcd100to250);
684     sample QCD250to500("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-250To500_7TeV-madgraph.root",((lumi/totEventsQCD250to500)*QCD250to500CrossSection),qcd250to500);
685     sample QCD500to1000("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-500To1000_7TeV-madgraph.root",((lumi/totEventsQCD500to1000)*QCD500to1000CrossSection),qcd500to1000);
686     sample QCD1000toInf("/scratch/fronga/MC_v1.29/QCD_TuneD6T_HT-1000_7TeV-madgraph.root",((lumi/totEventsQCD1000toInf)*QCD1000toInfCrossSection),qcd1000toInf);
687     sample LM4("/scratch/fronga/MC_v1.29/LM4_SUSY_sftsht_7TeV-pythia6.root",((lumi/(242322.0))*LM4CrossSection),lm4);
688    
689    
690    
691    
692    
693     // #------------
694    
695     cout<<"now I am armed ..."<<endl;
696     // #--- initialized histos
697    
698     createHisto(data);
699     createHisto(allB);
700     createHisto(zjets);
701     createHisto(ttbar);
702     createHisto(zjets_un);
703     createHisto(lm4);
704     createHisto(wjets);
705     createHisto(znunu);
706     createHisto(vvjets);
707     createHisto(singletops);
708     createHisto(singletopt);
709     createHisto(singletopu);
710     createHisto(qcd100to250);
711     createHisto(qcd250to500);
712     createHisto(qcd500to1000);
713     createHisto(qcd1000toInf);
714     createHisto(BnoZ);
715     createHisto(BwithZ);
716     createHisto(BandS);
717    
718     vector<sample> samples; // add samples to be projected (CAPS)
719     samples.push_back(Data);
720     samples.push_back(Zjets);
721     samples.push_back(TTbar);
722     samples.push_back(Zjets_un);
723     samples.push_back(LM4);
724     samples.push_back(WJets);
725     samples.push_back(Znunu);
726     samples.push_back(VVJets);
727     samples.push_back(SingleTopS);
728     samples.push_back(SingleTopT);
729     samples.push_back(SingleTopU);
730     samples.push_back(QCD100to250);
731     samples.push_back(QCD250to500);
732     samples.push_back(QCD500to1000);
733     samples.push_back(QCD1000toInf);
734    
735     // #--- define target sample build up from source samples
736     allBList.push_back(ttbar);
737     allBList.push_back(singletops);
738     allBList.push_back(singletopt);
739     allBList.push_back(singletopu);
740     allBList.push_back(znunu);
741     allBList.push_back(qcd100to250);
742     allBList.push_back(qcd250to500);
743     allBList.push_back(qcd500to1000);
744     allBList.push_back(qcd1000toInf);
745     allBList.push_back(wjets);
746     allBList.push_back(vvjets);
747     allBList.push_back(zjets);
748    
749     // #--- Start analysis
750    
751    
752     float lowlimit=2.5;
753     float highlimit=2.5;
754     bool DoNotPlot=false;
755     /*
756     bool electronChannel = false;
757     string numberOfJets = "ij0";
758     string legendHeader="";
759     if(electronChannel)
760     {
761     fileNamePreffix = "V1_zee_"+numberOfJets;
762     legendHeader = "Z(e^{+}e^{-})";
763     if(numberOfJets=="ij2")legendHeader=legendHeader+", nJets == 1";
764    
765     }else
766     {
767     fileNamePreffix = "V1_zmumu_"+numberOfJets;
768     legendHeader = "Z(#mu^{+}#mu^{-})";
769     if(numberOfJets=="ij2")legendHeader=legendHeader+", nJets == 1 ";
770     }
771     */
772     // TCut basiccut("(passed_triggers||genZPt>0)&&abs(pfJetEta[0])<1.2&&abs(pfJetEta[1])<1.2&&pfJetPt[0]>30&&pfJetPt[1]>30");
773     // TCut basiccut("pfJetPt[0]/pfJetScale[0]>30 && pfJetPt[1]/pfJetScale[1]>30 && abs(pfJetEta[0])<1.2&&abs(pfJetEta[1])<1.2");
774     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])");// Step 0
775     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||genZPt>0)");// Step 1
776     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||genZPt>0)&&abs(pfJetEta[0])<1.2&&abs(pfJetEta[1])<1.2"); // Step 2
777     TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||genZPt>0)"); // Step 2
778     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||genZPt>0)&&abs(pfJetEta[0])<1.2&&abs(pfJetEta[1])<1.2&&(Max$(abs(pfJetEta))<3.0)"); //Step 3
779     // TCut basiccut("(pfJetGoodNum>=2&&pfJetGoodID[0])&&(pfJetGoodNum>=2&&pfJetGoodID[1])&&(passed_triggers||genZPt>0)&&(Max$(abs(pfJetEta))<3.0)"); //Step 4
780     // TCut basiccut("mll>-1"); //if you don't want to require triggers, you can activate this basic cut :-)
781     TCut DiElectron("(id1+1)*(id2+1)*ch1*ch2==-1 && abs(eta1)<2.4 && abs(eta2)<2.4"&&basiccut);
782     TCut DiMuon("(id1+1)*(id2+1)*ch1*ch2==-4 && abs(eta1)<2.4 && abs(eta2)<2.4"&&basiccut);
783     // TCut OSSF("(id1+1)*(id2+1)*ch1*ch2!=-2 && abs(eta1)<2.4 && abs(eta2)<2.4&&id1==0");
784     TCut OSSF("(id1+1)*(id2+1)*ch1*ch2!=-2 && abs(eta1)<2.4 && abs(eta2)<2.4"&&basiccut);
785     TCut OSOF("(id1+1)*(id2+1)*ch1*ch2==-2 && abs(eta1)<2.4 && abs(eta2)<2.4"&&basiccut);
786     TCut DiLeptonMass("abs(mll-91.2)<20&&(passed_triggers||!passed_triggers||genZPt>0)");
787     // TCut DiLeptonMass("mll>10");
788     TCut Sideband("(mll>50 && mll<60) || (mll>120 && mll<160)"&&basiccut);
789     TCut METCut((METtype+">=60").c_str());
790     TCut ZPtCut("");
791     // TCut NJets("pfJetGoodNum>=2");
792     TCut NJets("(pfJetGoodNum>=3)"&&basiccut);
793     // TCut NJets("(pfJetGoodNum>=2)"&&basiccut);
794     // TCut NJets("jetNum>=2 && pt>40");
795     // TCut NJets("jetNum>=2 && pt>40");
796     // TCut NJets("pfJetGoodNum==1");
797    
798     string numberOfJets = "";
799     string legendHeader="";
800     legendHeader = "";
801     // float osofweight=0.5;// set this to 0.5 if you only consider ee (or mm) and not both ee & mm !!!
802     float osofweight=1.0;
803     // legendHeader = "m(l^{+}l^{-})>10 GeV && nJets #geq 2";
804     fileNamePreffix = "V3_2011_ll_incl";
805     TPaveText *eventSelectionPaveText = new TPaveText(0.27, 0.95, 0.77, 1.0,"blNDC");
806     eventSelectionPaveText->SetFillStyle(4000);
807     eventSelectionPaveText->SetFillColor(kWhite);
808     eventSelectionPaveText->SetTextFont(42);
809     // eventSelectionPaveText->SetBorderSize(1);
810     // eventSelectionPaveText->SetTextColor(kMagenta+2);
811     eventSelectionPaveText->AddText((legendHeader).c_str());
812    
813    
814    
815    
816    
817    
818    
819     // #--- construct the JZB distribution
820     for(int i=0;i<int(samples.size());i++) // project each sample in the proper histogram
821     {
822     samples[i].tree_->Project(jzb_ossf[samples[i].index_]->GetName(),JZBtype.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
823     samples[i].tree_->Project(jzb_ossfee[samples[i].index_]->GetName(),JZBtype.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && "id1==0")*samples[i].weight_*"weight");
824     samples[i].tree_->Project(jzb_ossfmm[samples[i].index_]->GetName(),JZBtype.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && "id1==1")*samples[i].weight_*"weight");
825     samples[i].tree_->Project(jzb_osof[samples[i].index_]->GetName(),JZBtype.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
826     jzb_osof[samples[i].index_]->Scale(osofweight);
827     }
828     // Add all MC
829     jzb_ossf[allB]= mc_th1f(jzb_ossf);
830     jzb_ossfee[allB]= mc_th1f(jzb_ossfee);
831     jzb_ossfmm[allB]= mc_th1f(jzb_ossfmm);
832     jzb_osof[allB]= mc_th1f(jzb_osof);
833    
834     // #--- apply subtraction
835     float subtractionWeight = -1.0;
836     float elefficiency=0.5170402668;
837     float muefficiency=0.4829597332;
838     jzb_subtracted[allB] = (TH1F*)jzb_ossf[allB]->Clone();
839     jzb_subtracted[allB]->Add(jzb_osof[allB],subtractionWeight);
840     jzb_subtractedee[allB] = (TH1F*)jzb_ossf[allB]->Clone();
841     jzb_subtractedee[allB]->Add(jzb_osof[allB],elefficiency*subtractionWeight);
842     jzb_subtractedmm[allB] = (TH1F*)jzb_ossf[allB]->Clone();
843     jzb_subtractedmm[allB]->Add(jzb_osof[allB],muefficiency*subtractionWeight);
844     for(int i=0;i<int(samples.size());i++) // project each sample in the proper histogram
845     {
846     jzb_subtracted[samples[i].index_] = (TH1F*)jzb_ossf[samples[i].index_]->Clone();
847     jzb_subtracted[samples[i].index_]->Add(jzb_osof[samples[i].index_], subtractionWeight);
848     jzb_subtractedee[samples[i].index_] = (TH1F*)jzb_ossf[samples[i].index_]->Clone();
849     jzb_subtractedee[samples[i].index_]->Add(jzb_osof[samples[i].index_], elefficiency*subtractionWeight);
850     jzb_subtractedmm[samples[i].index_] = (TH1F*)jzb_ossf[samples[i].index_]->Clone();
851     jzb_subtractedmm[samples[i].index_]->Add(jzb_osof[samples[i].index_], muefficiency*subtractionWeight);
852     }
853     jzb_subtracted[BandS] = (TH1F*)jzb_ossf[BandS]->Clone();
854     jzb_subtracted[BandS]->Add(jzb_osof[BandS], subtractionWeight);
855     jzb_subtractedee[BandS] = (TH1F*)jzb_ossf[BandS]->Clone();
856     jzb_subtractedee[BandS]->Add(jzb_osof[BandS], elefficiency*subtractionWeight);
857     jzb_subtractedmm[BandS] = (TH1F*)jzb_ossf[BandS]->Clone();
858     jzb_subtractedmm[BandS]->Add(jzb_osof[BandS], muefficiency*subtractionWeight);
859    
860     // #--- find JZB zero in MC
861     f1MC = new TF1("f1MC","gaus",-40,40);
862     float peakMCError=0.;
863     float MCsigma=0;
864     float peakMC = get_peak(jzb_subtracted[allB],peakMCError,MCsigma,f1MC, lowlimit, highlimit);
865     cout << "peakMC = " << peakMC << " +- " <<peakMCError<< " chi2 = " << f1MC->GetChisquare() << " NDF = " << f1MC->GetNDF() << " chi2/NDF = " << f1MC->GetChisquare()/f1MC->GetNDF() << endl;
866    
867     // #--- find JZB zero in Data
868     f1Data = new TF1("f1Data","gaus",-40,40);
869     float peakDataError=0.;
870     float datasigma=0;
871     float peakData = get_peak(jzb_subtracted[data],peakDataError,datasigma,f1Data, lowlimit, highlimit);
872     cout << "peakData = " << peakData << " +- " <<peakDataError<< " chi2 = " << f1Data->GetChisquare() << " NDF = " << f1Data->GetNDF() << " chi2/NDF = " << f1Data->GetChisquare()/f1Data->GetNDF() << endl;
873     /* cout << "ARTIFICIALLY OVERRIDDEN MU!" << endl;
874     float peakData=-1.3;
875     float peakDataError=0.6;*/
876    
877    
878     // #--- create new JZB binning
879     float ij2zeroPoint=0.;
880     float ij2zeroPointMisCalibP=0.;
881     float ij2zeroPointMisCalibN=0.;
882    
883    
884     float ij2zeroPointMC=peakMC;
885     float ij2zeroPointMisCalibPMC=peakMC+peakMCError;
886     float ij2zeroPointMisCalibNMC=peakMC-peakMCError;
887    
888    
889     float ij2zeroPointData=peakData;
890     float ij2zeroPointMisCalibPData=peakData+peakDataError;
891     float ij2zeroPointMisCalibNData=peakData-peakDataError;
892    
893    
894    
895    
896    
897     for(int i=0;i<int(samples.size());i++) // project each sample in the proper histogram
898     {
899    
900     if(i==data) // use the peak from the Data in the Data
901     {
902     ij2zeroPoint = ij2zeroPointData;
903     ij2zeroPointMisCalibP = ij2zeroPointMisCalibPData;
904     ij2zeroPointMisCalibN = ij2zeroPointMisCalibNData;
905     }else
906     {
907     ij2zeroPoint = ij2zeroPointMC;
908     ij2zeroPointMisCalibP = ij2zeroPointMisCalibPMC;
909     ij2zeroPointMisCalibN = ij2zeroPointMisCalibNMC;
910     }
911    
912     string JZBprime = JZBtype+"-"+any2string(ij2zeroPoint);
913     string sJZBprime = "(GausRandom("+JZBtype+"-"+any2string(ij2zeroPoint)+",5))";
914     string absJZBprime = "abs("+JZBtype+"-"+any2string(ij2zeroPoint)+")";
915     TCut JZBPrimeP((JZBtype+"-"+any2string(ij2zeroPoint)+">0").c_str()) ;
916     TCut JZBPrimeN((JZBtype+"-"+any2string(ij2zeroPoint)+"<0").c_str()) ;
917    
918    
919     string JZBprimeMisCalibP = JZBtype+"-"+any2string(ij2zeroPointMisCalibP);
920     string absJZBprimeMisCalibP = "abs("+JZBtype+"-"+any2string(ij2zeroPointMisCalibP)+")";
921     TCut JZBPrimeMisCalibPP((JZBtype+"-"+any2string(ij2zeroPointMisCalibP)+">0").c_str()) ;
922     TCut JZBPrimeMisCalibPN((JZBtype+"-"+any2string(ij2zeroPointMisCalibP)+"<0").c_str()) ;
923    
924     string JZBprimeMisCalibN = JZBtype+"-"+any2string(ij2zeroPointMisCalibN);
925     string absJZBprimeMisCalibN = "abs("+JZBtype+"-"+any2string(ij2zeroPointMisCalibN)+")";
926     TCut JZBPrimeMisCalibNP((JZBtype+"-"+any2string(ij2zeroPointMisCalibN)+">0").c_str()) ;
927     TCut JZBPrimeMisCalibNN((JZBtype+"-"+any2string(ij2zeroPointMisCalibN)+"<0").c_str()) ;
928    
929     if(i==data || i==allB)cout<<"sample: "<<i<<" ij2zeroPoint ="<<ij2zeroPoint<<endl;
930     if(i==data || i==allB)cout<<"sample: "<<i<<" ij2zeroPointMisCalibP ="<<ij2zeroPointMisCalibP<<endl;
931     if(i==data || i==allB)cout<<"sample: "<<i<<" ij2zeroPointMisCalibN ="<<ij2zeroPointMisCalibN<<endl;
932    
933     if(!DoNotPlot)
934     {
935     // #--- inclusive 2 lepton selection
936     samples[i].tree_->Project(mll[samples[i].index_]->GetName(),"mll",(OSSF)*samples[i].weight_*"weight");
937    
938     // #--- inclusive 2 lepton selection + Z mass requirement
939     samples[i].tree_->Project(J1Pt_ossf[samples[i].index_]->GetName(),"jetpt[0]",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
940     samples[i].tree_->Project(J1Eta_ossf[samples[i].index_]->GetName(),"jeteta[0]",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
941     samples[i].tree_->Project(ZPt_ossf_incl[samples[i].index_]->GetName(),"pt",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
942     samples[i].tree_->Project(leptonPt[samples[i].index_]->GetName(),"pt1",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
943     samples[i].tree_->Project(("+"+string(leptonPt[samples[i].index_]->GetName())).c_str(),"pt2",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
944     samples[i].tree_->Project(leptonEta[samples[i].index_]->GetName(),"eta1",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
945     samples[i].tree_->Project(("+"+string(leptonEta[samples[i].index_]->GetName())).c_str(),"eta2",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
946     samples[i].tree_->Project(nJets[samples[i].index_]->GetName(),"pfJetGoodNum",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
947     samples[i].tree_->Project(nJets285[samples[i].index_]->GetName(),"pfJetGoodNum285",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
948     samples[i].tree_->Project(nJets315[samples[i].index_]->GetName(),"pfJetGoodNum315",(OSSF && DiLeptonMass)*samples[i].weight_*"weight");
949    
950     // #--- jet selection + 2 lepton selection + Z mass requirement
951     samples[i].tree_->Project(JPt_ossf[samples[i].index_]->GetName(),"sumJetPt[1]",(OSSF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
952     samples[i].tree_->Project(ZPt_ossf[samples[i].index_]->GetName(),"pt",(OSSF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
953     }
954    
955     // #--- jet selection + 2 lepton selection + Z mass requirement ; plots
956     samples[i].tree_->Project(mll_ossf[samples[i].index_]->GetName(),"mll",(OSSF && NJets && ZPtCut)*samples[i].weight_*"weight");
957     samples[i].tree_->Project(mll_osof[samples[i].index_]->GetName(),"mll",(OSOF && NJets && ZPtCut)*samples[i].weight_*"weight");
958     mll_osof[samples[i].index_]->Scale(osofweight);
959    
960     samples[i].tree_->Project(mll_ee[samples[i].index_]->GetName(),"mll",(DiElectron && NJets && ZPtCut)*samples[i].weight_*"weight");
961     samples[i].tree_->Project(mll_mumu[samples[i].index_]->GetName(),"mll",(DiMuon && NJets && ZPtCut)*samples[i].weight_*"weight");
962    
963     // #--- jzb distributions
964     samples[i].tree_->Project(jzbprime_ossf[samples[i].index_]->GetName(),JZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
965     samples[i].tree_->Project(sjzbprime_ossf[samples[i].index_]->GetName(),sJZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
966     samples[i].tree_->Project(jzbprime_osof[samples[i].index_]->GetName(),JZBprime.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut)*samples[i].weight_*"weight");
967     jzbprime_osof[samples[i].index_]->Scale(osofweight);
968    
969     samples[i].tree_->Project(jzbprime_ossf_JZBP[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeP)*samples[i].weight_*"weight");
970     samples[i].tree_->Project(jzbprime_ossf_JZBPfin[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeP)*samples[i].weight_*"weight");
971     samples[i].tree_->Project(jzbprime_ossf_JZBN[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeN)*samples[i].weight_*"weight");
972     samples[i].tree_->Project(jzbprime_ossf_JZBNfin[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeN)*samples[i].weight_*"weight");
973     samples[i].tree_->Project(jzbprime_osof_JZBP[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeP)*samples[i].weight_*"weight");
974     samples[i].tree_->Project(jzbprime_osof_JZBPfin[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeP)*samples[i].weight_*"weight");
975     jzbprime_osof_JZBP[samples[i].index_]->Scale(osofweight);
976     jzbprime_osof_JZBPfin[samples[i].index_]->Scale(osofweight);
977     samples[i].tree_->Project(jzbprime_osof_JZBN[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeN)*samples[i].weight_*"weight");
978     samples[i].tree_->Project(jzbprime_osof_JZBNfin[samples[i].index_]->GetName(),absJZBprime.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeN)*samples[i].weight_*"weight");
979     jzbprime_osof_JZBN[samples[i].index_]->Scale(osofweight);
980     jzbprime_osof_JZBNfin[samples[i].index_]->Scale(osofweight);
981    
982     samples[i].tree_->Project(jzbprime_ossf_JZBP_MisCalibP[samples[i].index_]->GetName(),absJZBprimeMisCalibP.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibPP)*samples[i].weight_*"weight");
983     samples[i].tree_->Project(jzbprime_ossf_JZBN_MisCalibP[samples[i].index_]->GetName(),absJZBprimeMisCalibP.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibPN)*samples[i].weight_*"weight");
984     samples[i].tree_->Project(jzbprime_osof_JZBP_MisCalibP[samples[i].index_]->GetName(),absJZBprimeMisCalibP.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibPP)*samples[i].weight_*"weight");
985     jzbprime_osof_JZBP_MisCalibP[samples[i].index_]->Scale(osofweight);
986     samples[i].tree_->Project(jzbprime_osof_JZBN_MisCalibP[samples[i].index_]->GetName(),absJZBprimeMisCalibP.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibPN)*samples[i].weight_*"weight");
987     jzbprime_osof_JZBN_MisCalibP[samples[i].index_]->Scale(osofweight);
988    
989     samples[i].tree_->Project(jzbprime_ossf_JZBP_MisCalibN[samples[i].index_]->GetName(),absJZBprimeMisCalibN.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibNP)*samples[i].weight_*"weight");
990     samples[i].tree_->Project(jzbprime_ossf_JZBN_MisCalibN[samples[i].index_]->GetName(),absJZBprimeMisCalibN.c_str(),(OSSF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibNN)*samples[i].weight_*"weight");
991     samples[i].tree_->Project(jzbprime_osof_JZBP_MisCalibN[samples[i].index_]->GetName(),absJZBprimeMisCalibN.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibNP)*samples[i].weight_*"weight");
992     jzbprime_osof_JZBP_MisCalibN[samples[i].index_]->Scale(osofweight);
993     samples[i].tree_->Project(jzbprime_osof_JZBN_MisCalibN[samples[i].index_]->GetName(),absJZBprimeMisCalibN.c_str(),(OSOF && DiLeptonMass && NJets && ZPtCut && JZBPrimeMisCalibNN)*samples[i].weight_*"weight");
994     jzbprime_osof_JZBN_MisCalibN[samples[i].index_]->Scale(osofweight);
995     }
996    
997     // #--- sum up histos
998     for(int i=0;i<int(allTH1F.size());i++)
999     {
1000     if(jzb_ossf[allB]==allTH1F[i][allB])continue; // pass-by this sample (done already)
1001     if(jzb_osof[allB]==allTH1F[i][allB])continue; // pass-by this sample (done already)
1002    
1003     allTH1F[i][allB] = mc_th1f(allTH1F[i]);
1004     }
1005    
1006    
1007     cout<<"samples have been digested ..."<<endl;
1008    
1009    
1010    
1011     gStyle->SetMarkerSize(1.2);
1012    
1013    
1014     // #--- start plotting
1015     canvasCounter =0;
1016    
1017    
1018     string overallTitle="CMS preliminary 2011";
1019     stringstream Cmenstream;
1020     Cmenstream << "#sqrt{s}=7TeV, L = " << lumi << " pb^{-1}";
1021     // string CMen="#sqrt{s}=7TeV, L = 34 pb^{-1}";
1022     string CMen=Cmenstream.str();
1023     paveTitle = new TPaveText(0.6,0.8,0.95,0.93,"blNDC");
1024     paveTitle->SetFillStyle(4000);
1025     paveTitle->SetFillColor(kWhite);
1026     //paveTitle->SetBorderSize(0.1);
1027     paveTitle->SetTextFont(42);
1028     paveTitle->AddText(overallTitle.c_str());
1029     paveTitle->AddText(CMen.c_str());
1030    
1031     paveTitleMCZjets = new TPaveText(0.6,0.8,0.95,0.93,"blNDC");
1032     paveTitleMCZjets->SetFillStyle(4000);
1033     paveTitleMCZjets->SetFillColor(kWhite);
1034     //paveTitleMCZjets->SetBorderSize(0.1);
1035     paveTitleMCZjets->SetTextFont(42);
1036     paveTitleMCZjets->AddText("Madgraph Z+jets mock data");
1037     paveTitleMCZjets->AddText("#sqrt{s}=7TeV, L = 900 pb^{-1}");
1038    
1039     paveTitleMC = new TPaveText(0.6,0.8,0.95,0.93,"blNDC");
1040     paveTitleMC->SetFillStyle(4000);
1041     paveTitleMC->SetFillColor(kWhite);
1042     //paveTitleMC->SetBorderSize(0.1);
1043     paveTitleMC->SetTextFont(42);
1044     paveTitleMC->AddText("CMS MC simulation");
1045     stringstream paveTitleMCtext;
1046     paveTitleMCtext << "#sqrt{s}=7TeV, L = "<<lumi<<" pb^{-1}";
1047     // paveTitleMC->AddText("#sqrt{s}=7TeV, L = 34 pb^{-1}");
1048     paveTitleMC->AddText(paveTitleMCtext.str().c_str());
1049    
1050    
1051     globalMCLeg = new TLegend(0.66,0.64,0.93,0.78);
1052     globalMCLeg->SetFillColor(kWhite);
1053     globalMCLeg->SetTextFont(42);
1054     // globalMCLeg->SetHeader(legendHeader.c_str());
1055     jzb_ossf[zjets]->SetFillColor(kYellow);
1056     jzb_ossf[ttbar]->SetFillColor(kMagenta+2);
1057     jzb_ossf[wjets]->SetFillColor(kBlack);
1058    
1059    
1060     jzb_ossf[ttbar]->SetFillColor(kMagenta+2);
1061     jzb_ossf[zjets]->SetFillColor(kYellow);
1062     jzb_ossf[wjets]->SetFillColor(kGray);
1063     jzb_ossf[znunu]->SetFillColor(kYellow);
1064     jzb_ossf[vvjets]->SetFillColor(kGreen);
1065     jzb_ossf[singletops]->SetFillColor(kBlue);
1066     jzb_ossf[singletopt]->SetFillColor(kBlue);
1067     jzb_ossf[singletopu]->SetFillColor(kBlue);
1068     jzb_ossf[qcd100to250]->SetFillColor(kPink);
1069    
1070    
1071     globalMCLeg->AddEntry(jzb_ossf[zjets],"Zjets","f");
1072     globalMCLeg->AddEntry(jzb_ossf[ttbar],"t#bar{t}","f");
1073     globalMCLeg->AddEntry(jzb_ossf[wjets],"WJets","f");
1074     globalMCLeg->AddEntry(jzb_ossf[vvjets],"DiBosons","f");
1075     globalMCLeg->AddEntry(jzb_ossf[singletops],"SingleTop","f");
1076     globalMCLeg->AddEntry(jzb_ossf[qcd100to250],"QCD","f");
1077    
1078    
1079     globalDataMCLeg = new TLegend(0.66,0.64,0.93,0.78);
1080     globalDataMCLeg->SetFillColor(kWhite);
1081     globalDataMCLeg->SetTextFont(42);
1082     globalDataMCLeg->AddEntry(jzb_ossf[zjets],"Zjets","f");
1083     globalDataMCLeg->AddEntry(jzb_ossf[ttbar],"t#bar{t}","f");
1084     globalDataMCLeg->AddEntry(jzb_ossf[wjets],"WJets","f");
1085     globalDataMCLeg->AddEntry(jzb_ossf[vvjets],"DiBosons","f");
1086     globalDataMCLeg->AddEntry(jzb_ossf[singletops],"SingleTop","f");
1087     globalDataMCLeg->AddEntry(jzb_ossf[qcd100to250],"QCD","f");
1088    
1089    
1090     // #--- always plot this!
1091     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_subtracted_mc"+any2string(canvasCounter)).c_str());
1092     can[canvasCounter]->cd();
1093     canvasCounter++;
1094     hs[canvasCounter-1] = mc_stack(jzb_subtracted);
1095     jzb_subtracted[allB]->SetMinimum(0);
1096     // setYMax(jzb_subtracted[allB],1.2);
1097     jzb_subtracted[allB]->GetXaxis()->SetRangeUser(-50,50);
1098     jzb_subtracted[allB]->SetMaximum(200);
1099     jzb_subtracted[allB]->SetMinimum(0);
1100     jzb_subtracted[allB]->Draw("hist");
1101     hs[canvasCounter-1]->Draw("hist,same");
1102     f1MC->Draw("same");
1103     paveTitle->Draw("same");
1104     jzb_subtracted[allB]->Draw("axis,same");
1105     TPaveText * paveFitMC = new TPaveText(0.16,0.67,0.40,0.76,"blNDC");
1106     paveFitMC->SetFillStyle(4000);
1107     paveFitMC->SetFillColor(kWhite);
1108     paveFitMC->SetTextFont(42);
1109     ostringstream firstLineMC;
1110     ostringstream secondLineMC;
1111     firstLineMC<< setprecision(2) << "#mu=" << f1MC->GetParameter(1) << "#pm" << f1MC->GetParError(1) << " (GeV)";
1112     secondLineMC<< setprecision(2) << "#sigma=" << f1MC->GetParameter(2) << "#pm" << f1MC->GetParError(2) << " (GeV)";
1113     paveFitMC->AddText(firstLineMC.str().c_str());
1114     paveFitMC->AddText(secondLineMC.str().c_str());
1115     paveFitMC->Draw("same");
1116     globalMCLeg->Draw("same");
1117     eventSelectionPaveText->Draw("same");
1118     can[canvasCounter-1]->Update();
1119     saveCanvas(canvasCounter-1,"jzb_subtracted_mc");
1120    
1121     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_subtracted_data"+any2string(canvasCounter)).c_str());
1122     can[canvasCounter]->cd();
1123     canvasCounter++;
1124     jzb_subtracted[data]->GetXaxis()->SetRangeUser(-50,50);
1125     jzb_subtracted[data]->SetMinimum(0);
1126     jzb_subtracted[data]->Draw("e1");
1127     f1Data->Draw("same");
1128     paveTitle->Draw("same");
1129     jzb_subtracted[data]->Draw("axis, same");
1130     TPaveText * paveFitData = new TPaveText(0.16,0.67,0.40,0.76,"blNDC");
1131     paveFitData->SetFillStyle(4000);
1132     paveFitData->SetFillColor(kWhite);
1133     paveFitData->SetTextFont(42);
1134     ostringstream firstLineData;
1135     ostringstream secondLineData;
1136     firstLineData<< setprecision(2) << "#mu=" << f1Data->GetParameter(1) << "#pm" << f1Data->GetParError(1) << " (GeV)";
1137     secondLineData<< setprecision(2) << "#sigma=" << f1Data->GetParameter(2) << "#pm" << f1Data->GetParError(2) << " (GeV)";
1138     paveFitData->AddText(firstLineData.str().c_str());
1139     paveFitData->AddText(secondLineData.str().c_str());
1140     paveFitData->Draw("same");
1141     leg[canvasCounter] = new TLegend(0.69,0.69,0.87,0.75);
1142     leg[canvasCounter] ->SetFillStyle(4000);
1143     leg[canvasCounter] ->SetFillColor(kWhite);
1144     leg[canvasCounter] ->SetTextFont(42);
1145     leg[canvasCounter] ->AddEntry(jzb_subtracted[data],"data","p");
1146     leg[canvasCounter] ->Draw("same");
1147     eventSelectionPaveText->Draw("same");
1148     can[canvasCounter-1]->Update();
1149     saveCanvas(canvasCounter-1,"jzb_subtracted_data");
1150    
1151    
1152     if(!DoNotPlot)
1153     {
1154    
1155     // #---
1156     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("J1Pt_ossf_c"+any2string(canvasCounter)).c_str());
1157     can[canvasCounter]->cd();
1158     can[canvasCounter]->SetLogy(1);
1159     canvasCounter++;
1160     J1Pt_ossf[data]->Draw("e1");
1161     hs[canvasCounter-1] = mc_stack(J1Pt_ossf);
1162     hs[canvasCounter-1]->Draw("same,hist");
1163     J1Pt_ossf[data]->Draw("e1,same");
1164     J1Pt_ossf[data]->Draw("axis,same");
1165     paveTitle->Draw("same");
1166     globalDataMCLeg->Draw("same");
1167     eventSelectionPaveText->Draw("same");
1168     can[canvasCounter-1]->Update();
1169     saveCanvas(canvasCounter-1);
1170    
1171     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("J1Eta_ossf_c"+any2string(canvasCounter)).c_str());
1172     can[canvasCounter]->cd();
1173     canvasCounter++;
1174     setYMax(J1Eta_ossf[data],1.4);
1175     J1Eta_ossf[data]->Draw("e1");
1176     hs[canvasCounter-1] = mc_stack(J1Eta_ossf);
1177     hs[canvasCounter-1]->Draw("same,hist");
1178     J1Eta_ossf[data]->Draw("e1,same");
1179     J1Eta_ossf[data]->Draw("axis,same");
1180     paveTitle->Draw("same");
1181     globalDataMCLeg->Draw("same");
1182     eventSelectionPaveText->Draw("same");
1183     can[canvasCounter-1]->Update();
1184     saveCanvas(canvasCounter-1);
1185    
1186     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("ZPt_ossf_incl_c"+any2string(canvasCounter)).c_str());
1187     can[canvasCounter]->cd();
1188     can[canvasCounter]->SetLogy(1);
1189     canvasCounter++;
1190     ZPt_ossf_incl[data]->Draw("e1");
1191     ZPt_ossf_incl[data]->SetMinimum(1);
1192     hs[canvasCounter-1] = mc_stack(ZPt_ossf_incl);
1193     hs[canvasCounter-1]->Draw("same,hist");
1194     ZPt_ossf_incl[data]->Draw("e1,same");
1195     ZPt_ossf_incl[data]->Draw("axis,same");
1196     paveTitle->Draw("same");
1197     globalDataMCLeg->Draw("same");
1198     eventSelectionPaveText->Draw("same");
1199     can[canvasCounter-1]->Update();
1200     saveCanvas(canvasCounter-1);
1201    
1202    
1203     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("leptonPt_c"+any2string(canvasCounter)).c_str());
1204     can[canvasCounter]->cd();
1205     canvasCounter++;
1206     setYMax(leptonPt[data],1.4);
1207     leptonPt[data]->GetXaxis()->SetRangeUser(0,100);
1208     leptonPt[data]->Draw("e1");
1209     hs[canvasCounter-1] = mc_stack(leptonPt);
1210     hs[canvasCounter-1]->Draw("same,hist");
1211     leptonPt[data]->Draw("e1,same");
1212     leptonPt[data]->Draw("axis,same");
1213     paveTitle->Draw("same");
1214     globalDataMCLeg->Draw("same");
1215     eventSelectionPaveText->Draw("same");
1216     can[canvasCounter-1]->Update();
1217     saveCanvas(canvasCounter-1);
1218    
1219    
1220     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("leptonEta_c"+any2string(canvasCounter)).c_str());
1221     can[canvasCounter]->cd();
1222     canvasCounter++;
1223     setYMax(leptonEta[data],2.0);
1224     leptonEta[data]->Draw("e1");
1225     hs[canvasCounter-1] = mc_stack(leptonEta);
1226     hs[canvasCounter-1]->Draw("same,hist");
1227     leptonEta[data]->Draw("e1,same");
1228     leptonEta[data]->Draw("axis,same");
1229     paveTitle->Draw("same");
1230     globalDataMCLeg->Draw("same");
1231     eventSelectionPaveText->Draw("same");
1232     can[canvasCounter-1]->Update();
1233     saveCanvas(canvasCounter-1);
1234    
1235    
1236     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("mll_c"+any2string(canvasCounter)).c_str());
1237     can[canvasCounter]->cd();
1238     can[canvasCounter]->SetLogy(1);
1239     canvasCounter++;
1240     setYMax(mll[data],2.2);
1241     mll[data]->GetXaxis()->SetRangeUser(50,150);
1242     mll[data]->GetYaxis()->SetRangeUser(1,mll[data]->GetMaximum()*1.2);
1243     mll[data]->Draw("e1");
1244     mll[data]->Draw("axis,same");
1245     hs[canvasCounter-1] = mc_stack(mll);
1246     hs[canvasCounter-1]->Draw("same,hist");
1247     mll[data]->Draw("e1,same");
1248     mll[data]->Draw("axis,same");
1249     paveTitle->Draw("same");
1250     globalDataMCLeg->Draw("same");
1251     eventSelectionPaveText->Draw("same");
1252     can[canvasCounter-1]->Update();
1253     saveCanvas(canvasCounter-1);
1254     }//end of DoNotPlot
1255    
1256     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_ossf_c"+any2string(canvasCounter)).c_str());
1257     can[canvasCounter]->cd();
1258     can[canvasCounter]->SetLogy(1);
1259    
1260     canvasCounter++;
1261     setYMax(jzbprime_ossf[data],1.4);
1262     // jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-150,150);
1263     jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-120,350);
1264     jzbprime_ossf[data]->SetMinimum(0.1);
1265     jzbprime_ossf[data]->SetMaximum(1500);
1266     jzbprime_ossf[data]->Draw("e1");
1267     hs[canvasCounter-1] = mc_stack(jzbprime_ossf);
1268     hs[canvasCounter-1]->Draw("same,hist");
1269     jzbprime_ossf[data]->Draw("e1,same");
1270     jzbprime_ossf[data]->Draw("axis,same");
1271     paveTitle->Draw("same");
1272     globalDataMCLeg->Draw("same");
1273     eventSelectionPaveText->Draw("same");
1274     can[canvasCounter-1]->Update();
1275     saveCanvas(canvasCounter-1);
1276    
1277     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("nJets_c"+any2string(canvasCounter)).c_str());
1278     can[canvasCounter]->cd();
1279     can[canvasCounter]->SetLogy(1);
1280     canvasCounter++;
1281     nJets[data]->GetXaxis()->SetRangeUser(0,7);
1282     nJets[data]->Draw("e1");
1283     hs[canvasCounter-1] = mc_stack(nJets);
1284     hs[canvasCounter-1]->Draw("same,hist");
1285     nJets[data]->Draw("e1,same");
1286     nJets[data]->Draw("axis,same");
1287     paveTitle->Draw("same");
1288     globalDataMCLeg->Draw("same");
1289     eventSelectionPaveText->Draw("same");
1290     can[canvasCounter-1]->Update();
1291     saveCanvas(canvasCounter-1);
1292    
1293    
1294    
1295     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_ossf_c"+any2string(canvasCounter)).c_str());
1296     can[canvasCounter]->cd();
1297     can[canvasCounter]->SetLogy(1);
1298    
1299     canvasCounter++;
1300     setYMax(jzbprime_ossf[data],1.4);
1301     // jzbprime_ossf[data]->SetMinimum(0.0001);
1302     // jzbprime_ossf[data]->SetMaximum(1);
1303     jzbprime_ossf[data]->DrawNormalized("e1");
1304     hs[canvasCounter-1] = mc_normalized_stack(sjzbprime_ossf);
1305     hs[canvasCounter-1]->Draw("same,hist");
1306     jzbprime_ossf[data]->DrawNormalized("e1,same");
1307     jzbprime_ossf[data]->DrawNormalized("axis,same");
1308     paveTitle->Draw("same");
1309     globalDataMCLeg->Draw("same");
1310     eventSelectionPaveText->Draw("same");
1311     can[canvasCounter-1]->Update();
1312     saveCanvas(canvasCounter-1,"Normalized_SJZB");
1313    
1314    
1315     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_normalized_ossf_SHAPECOMPARISON_c"+any2string(canvasCounter)).c_str());
1316     can[canvasCounter]->cd();
1317     can[canvasCounter]->SetLogy(1);
1318    
1319     canvasCounter++;
1320     setYMax(jzbprime_ossf[data],1.4);
1321     // setYMax(jzbprime_ossf[data],1.4);
1322     // jzbprime_ossf[data]->SetMinimum(0.0001);
1323     // jzbprime_ossf[data]->SetMaximum(1);
1324     // jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-100,200);
1325     jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-120,350);
1326     // jzbprime_ossf[data]->DrawNormalized("e1");
1327     jzbprime_ossf[data]->SetMinimum(0.1);
1328     jzbprime_ossf[data]->SetMaximum(200);
1329     // jzbprime_ossf[data]->SetMaximum(1.2*jzbprime_ossf[data]->GetMaximum());
1330     jzbprime_ossf[data]->Draw("e1");
1331     // hs[canvasCounter-1] = mc_normalized_stack(jzbprime_ossf);
1332     hs[canvasCounter-1] = mc_normalizedtodata_stack(jzbprime_ossf);
1333    
1334     hs[canvasCounter-1]->Draw("same,hist");
1335     // jzbprime_ossf[data]->DrawNormalized("e1,same");
1336     // jzbprime_ossf[data]->DrawNormalized("axis,same");
1337     jzbprime_ossf[data]->Draw("e1,same");
1338     jzbprime_ossf[data]->Draw("axis,same");
1339     paveTitle->Draw("same");
1340     globalDataMCLeg->Draw("same");
1341     eventSelectionPaveText->Draw("same");
1342     can[canvasCounter-1]->Update();
1343     saveCanvas(canvasCounter-1,"NORMALIZED_OSSF_JZB__SHAPECOMPARISON");
1344    
1345     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_osof_normalized_SHAPECOMPARISON_c"+any2string(canvasCounter)).c_str());
1346     can[canvasCounter]->cd();
1347     can[canvasCounter]->SetLogy(1);
1348    
1349     canvasCounter++;
1350     // setYMax(jzbprime_osof[data],1.4);
1351     jzbprime_osof[data]->SetMaximum(40);
1352     // jzbprime_osof[data]->GetYaxis()->SetRangeUser(0.0001,1);
1353     // jzbprime_osof[data]->DrawNormalized("e1");
1354     setYMax(jzbprime_osof[data],1.4);
1355     jzbprime_osof[data]->SetMinimum(0.1);
1356     jzbprime_osof[data]->Draw("e1");
1357     // jzbprime_osof[data]->GetXaxis()->SetRangeUser(-100,200);
1358     jzbprime_osof[data]->GetXaxis()->SetRangeUser(-120,350);
1359     // hs[canvasCounter-1] = mc_normalized_stack(jzbprime_osof);
1360     hs[canvasCounter-1] = mc_normalizedtodata_stack(jzbprime_osof);
1361     hs[canvasCounter-1]->Draw("same,hist");
1362     // jzbprime_osof[data]->DrawNormalized("e1,same");
1363     // jzbprime_osof[data]->DrawNormalized("axis,same");
1364     jzbprime_osof[data]->Draw("e1,same");
1365     jzbprime_osof[data]->Draw("axis,same");
1366     paveTitle->Draw("same");
1367     globalDataMCLeg->Draw("same");
1368     eventSelectionPaveText->Draw("same");
1369     can[canvasCounter-1]->Update();
1370     saveCanvas(canvasCounter-1,"NORMALIZED_OSOF_JZB__SHAPECOMPARISON");
1371     /*
1372     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_osof_unnormalized_c"+any2string(canvasCounter)).c_str());
1373     can[canvasCounter]->cd();
1374     can[canvasCounter]->SetLogy(1);
1375    
1376     canvasCounter++;
1377     // setYMax(jzbprime_osof[data],1.4);
1378     // jzbprime_osof[data]->GetYaxis()->SetRangeUser(0.0001,1);
1379     // jzbprime_osof[data]->GetXaxis()->SetRangeUser(-100,200);
1380     jzbprime_osof[data]->GetXaxis()->SetRangeUser(-120,350);
1381     jzbprime_osof[data]->GetYaxis()->SetRangeUser(0.1,10);
1382     jzbprime_osof[data]->Draw("e1");
1383     hs[canvasCounter-1] = mc_stack(jzbprime_osof);
1384     hs[canvasCounter-1]->Draw("same,hist");
1385     jzbprime_osof[data]->Draw("e1,same");
1386     jzbprime_osof[data]->Draw("axis,same");
1387     paveTitle->Draw("same");
1388     globalDataMCLeg->Draw("same");
1389     eventSelectionPaveText->Draw("same");
1390     can[canvasCounter-1]->Update();
1391     saveCanvas(canvasCounter-1,"UNNORMALIZED_OSOF_JZB");
1392    
1393    
1394     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_ossf_unnormalized_c"+any2string(canvasCounter)).c_str());
1395     can[canvasCounter]->cd();
1396     can[canvasCounter]->SetLogy(1);
1397    
1398     canvasCounter++;
1399     // setYMax(jzbprime_osof[data],1.4);
1400     // jzbprime_osof[data]->GetYaxis()->SetRangeUser(0.0001,1);
1401     // jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-100,200);
1402     jzbprime_ossf[data]->GetXaxis()->SetRangeUser(-120,350);
1403     jzbprime_ossf[data]->GetYaxis()->SetRangeUser(0.1,jzbprime_ossf[data]->GetMaximum()*1.2);
1404     jzbprime_ossf[data]->Draw("e1");
1405     hs[canvasCounter-1] = mc_stack(jzbprime_ossf);
1406     hs[canvasCounter-1]->Draw("same,hist");
1407     jzbprime_ossf[data]->Draw("e1,same");
1408     jzbprime_ossf[data]->Draw("axis,same");
1409     paveTitle->Draw("same");
1410     globalDataMCLeg->Draw("same");
1411     eventSelectionPaveText->Draw("same");
1412     can[canvasCounter-1]->Update();
1413     saveCanvas(canvasCounter-1,"UNNORMALIZED_OSSF_JZB");
1414    
1415     */
1416     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzbprime_osof_unnormalized_rebinned_SHAPECOMPARISON_c"+any2string(canvasCounter)).c_str());
1417     can[canvasCounter]->cd();
1418     can[canvasCounter]->SetLogy(1);
1419     TH1F *jzbprime_osof_rebin[20];
1420     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
1421     {
1422     int sampleIndex = allBList[counter];
1423     jzbprime_osof_rebin[sampleIndex]=(TH1F*)jzbprime_osof[counter]->Clone();
1424     jzbprime_osof_rebin[sampleIndex]->Rebin(4);
1425     cout << "For sample " << sampleIndex << " we have an integral of " << jzbprime_osof_rebin[sampleIndex]->Integral() << endl;
1426     }
1427     jzbprime_osof_rebin[data]=(TH1F*)jzbprime_osof[data]->Clone();
1428     jzbprime_osof_rebin[data]->Rebin(4);
1429     jzbprime_osof_rebin[zjets]->SetName("jzbprime_osof_rebinned_zjets");
1430     canvasCounter++;
1431     // setYMax(jzbprime_ossf[data],1.4);
1432     jzbprime_osof_rebin[data]->SetMaximum(40);
1433     jzbprime_osof_rebin[data]->SetMinimum(0.1);
1434     // jzbprime_osof_rebin[data]->GetXaxis()->SetRangeUser(-100,200);
1435     jzbprime_osof_rebin[data]->GetXaxis()->SetRangeUser(-120,350);
1436     jzbprime_osof_rebin[data]->Draw("e1");
1437     // hs[canvasCounter-1] = mc_stack(jzbprime_osof_rebin);
1438     // hs[canvasCounter-1] = mc_stack(jzbprime_osof,4);
1439     hs[canvasCounter-1] = mc_normalizedtodata_stack(jzbprime_osof,4);
1440     hs[canvasCounter-1]->Draw("same,hist");
1441     jzbprime_osof_rebin[data]->Draw("e1,same");
1442     jzbprime_osof_rebin[data]->Draw("axis,same");
1443     paveTitle->Draw("same");
1444     globalDataMCLeg->Draw("same");
1445     eventSelectionPaveText->Draw("same");
1446     can[canvasCounter-1]->Update();
1447     saveCanvas(canvasCounter-1,"UNNORMALIZED_OSOF_REBINNED_JZB__SHAPECOMPARISON");
1448    
1449     if(!DoNotPlot) {
1450     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("nJets285_c"+any2string(canvasCounter)).c_str());
1451     can[canvasCounter]->cd();
1452     can[canvasCounter]->SetLogy(1);
1453     canvasCounter++;
1454     nJets285[data]->GetXaxis()->SetRangeUser(0,7);
1455     nJets285[data]->Draw("e1");
1456     hs[canvasCounter-1] = mc_stack(nJets285);
1457     hs[canvasCounter-1]->Draw("same,hist");
1458     nJets285[data]->Draw("e1,same");
1459     nJets285[data]->Draw("axis,same");
1460     paveTitle->Draw("same");
1461     globalDataMCLeg->Draw("same");
1462     eventSelectionPaveText->Draw("same");
1463     can[canvasCounter-1]->Update();
1464     saveCanvas(canvasCounter-1,"nJets285");
1465    
1466     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("nJets315_c"+any2string(canvasCounter)).c_str());
1467     can[canvasCounter]->cd();
1468     can[canvasCounter]->SetLogy(1);
1469     canvasCounter++;
1470     nJets315[data]->GetXaxis()->SetRangeUser(0,7);
1471     nJets315[data]->Draw("e1");
1472     hs[canvasCounter-1] = mc_stack(nJets315);
1473     hs[canvasCounter-1]->Draw("same,hist");
1474     nJets315[data]->Draw("e1,same");
1475     nJets315[data]->Draw("axis,same");
1476     paveTitle->Draw("same");
1477     globalDataMCLeg->Draw("same");
1478     eventSelectionPaveText->Draw("same");
1479     can[canvasCounter-1]->Update();
1480     saveCanvas(canvasCounter-1,"nJets315");
1481    
1482     }// #--- end of DoNotPlot scope
1483    
1484     //TF1 * logPol = new TF1("logPol","[0]*([1] + x)**[2]",0,1000);Bprediction->Fit(logPol,"LL","",25,500);
1485    
1486     // calculate e/mu efficiency ratio
1487     float mllmin;
1488     float mllmax;
1489     float Nee;
1490     float Nmumu;
1491     float eff_ratioSquare_allB;
1492     float eff_ratioSquare_data;
1493     float eff_ratioSquare_Error_allB;
1494     float eff_ratioSquare_Error_data;
1495     float eff_ratio_allB;
1496     float eff_ratio_data;
1497     float eff_ratio_Error_allB;
1498     float eff_ratio_Error_data;
1499     float emuNorm_allB;
1500     float emuNormError_allB;
1501     float emuNorm_data;
1502     float emuNormError_data;
1503    
1504     mllmin = 70; mllmax=110 ;
1505    
1506    
1507     Nee = histIntegral(mll_ee[allB],mllmin,mllmax);
1508     Nmumu = histIntegral(mll_mumu[allB],mllmin,mllmax);
1509     eff_ratioSquare_allB = Nee/Nmumu;
1510     eff_ratioSquare_Error_allB = computeRatioError(Nee,sqrt(Nee),Nmumu,sqrt(Nmumu));
1511     eff_ratio_allB = sqrt(eff_ratioSquare_allB);
1512     eff_ratio_Error_allB = 0.5*(1/sqrt(eff_ratioSquare_allB))*eff_ratioSquare_Error_allB;
1513     emuNorm_allB = 0.5*(eff_ratio_allB + (1/eff_ratio_allB) );
1514     emuNormError_allB = ( 1 - (1/(eff_ratio_allB*eff_ratio_allB)) )*eff_ratio_Error_allB;
1515    
1516     cout << "eff_ratio_allB = " << eff_ratio_allB << " +- " << eff_ratio_Error_allB <<endl;
1517     cout << "emuNorm_allB = " << emuNorm_allB << " +- " << emuNormError_allB << endl;
1518    
1519    
1520     Nee = histIntegral(mll_ee[data],mllmin,mllmax);
1521     Nmumu = histIntegral(mll_mumu[data],mllmin,mllmax);
1522     eff_ratioSquare_data = Nee/Nmumu;
1523     eff_ratioSquare_Error_data = computeRatioError(Nee,sqrt(Nee),Nmumu,sqrt(Nmumu));
1524     eff_ratio_data = sqrt(eff_ratioSquare_data);
1525     eff_ratio_Error_data = 0.5*(1/sqrt(eff_ratioSquare_data))*eff_ratioSquare_Error_data;
1526     emuNorm_data = 0.5*(eff_ratio_data + (1/eff_ratio_data) );
1527     emuNormError_data = ( 1 - (1/(eff_ratio_data*eff_ratio_data)) )*eff_ratio_Error_data;
1528    
1529     cout << "eff_ratio_data = " << eff_ratio_data << " +- " << eff_ratio_Error_data <<endl;
1530     cout << "emuNorm_data = " << emuNorm_data << " +- " << emuNormError_data << endl;
1531    
1532    
1533    
1534    
1535     //
1536    
1537     float JZBmaxMC = 300;
1538     float JZBmaxData = 200;
1539     float range=0;
1540     bool plotSys = false;
1541    
1542     float rMC=1.;
1543     float rMCSysP=1.00;
1544     float rMCSysN=1.00;
1545    
1546     Bpred[allB] = (TH1F*)jzbprime_ossf_JZBN[allB]->Clone();
1547     Bpred[allB]->Add(jzbprime_osof_JZBN[allB],-rMC);
1548     Bpred[allB]->Add(jzbprime_osof_JZBP[allB],+rMC);
1549    
1550     BpredSysP[allB] = (TH1F*)jzbprime_ossf_JZBN_MisCalibP[allB]->Clone();
1551     BpredSysP[allB]->Add(jzbprime_osof_JZBN_MisCalibP[allB],-rMCSysP);
1552     BpredSysP[allB]->Add(jzbprime_osof_JZBP_MisCalibP[allB],+rMCSysP);
1553    
1554     BpredSysN[allB] = (TH1F*)jzbprime_ossf_JZBN_MisCalibN[allB]->Clone();
1555     BpredSysN[allB]->Add(jzbprime_osof_JZBN_MisCalibN[allB],-rMCSysN);
1556     BpredSysN[allB]->Add(jzbprime_osof_JZBP_MisCalibN[allB],+rMCSysN);
1557    
1558     Bpredfin[allB] = (TH1F*)jzbprime_ossf_JZBNfin[allB]->Clone();
1559     Bpredfin[allB]->Add(jzbprime_osof_JZBNfin[allB],-rMC);
1560     Bpredfin[allB]->Add(jzbprime_osof_JZBPfin[allB],+rMC);
1561    
1562     // JZBratio[allB] = histRatio(jzbprime_ossf_JZBP[allB], Bpred[allB], allB);
1563     JZBratio[allB] = histRatiofin(jzbprime_ossf_JZBPfin[allB], Bpredfin[allB], allB);
1564    
1565     pol0[allB] = new TF1(("s"+any2string(allB)).c_str(),"[0]",0,1000);
1566     JZBratio[allB]->Fit(pol0[allB],"Q0R","",0,30);
1567    
1568     Bpred[BandS] = (TH1F*)jzbprime_ossf_JZBN[BandS]->Clone();
1569     Bpred[BandS]->Add(jzbprime_osof_JZBN[BandS],-rMC);
1570     Bpred[BandS]->Add(jzbprime_osof_JZBP[BandS],+rMC);
1571    
1572     Bpredfin[BandS] = (TH1F*)jzbprime_ossf_JZBNfin[BandS]->Clone();
1573     Bpredfin[BandS]->Add(jzbprime_osof_JZBNfin[BandS],-rMC);
1574     Bpredfin[BandS]->Add(jzbprime_osof_JZBPfin[BandS],+rMC);
1575    
1576     BpredSysP[BandS] = (TH1F*)jzbprime_ossf_JZBN_MisCalibP[BandS]->Clone();
1577     BpredSysP[BandS]->Add(jzbprime_osof_JZBN_MisCalibP[BandS],-rMCSysP);
1578     BpredSysP[BandS]->Add(jzbprime_osof_JZBP_MisCalibP[BandS],+rMCSysP);
1579    
1580     BpredSysN[BandS] = (TH1F*)jzbprime_ossf_JZBN_MisCalibN[BandS]->Clone();
1581     BpredSysN[BandS]->Add(jzbprime_osof_JZBN_MisCalibN[BandS],-rMCSysN);
1582     BpredSysN[BandS]->Add(jzbprime_osof_JZBP_MisCalibN[BandS],+rMCSysN);
1583    
1584    
1585     JZBratio[BandS] = histRatio(jzbprime_ossf_JZBPfin[BandS], Bpredfin[BandS], BandS);
1586     pol0[BandS] = new TF1(("s"+any2string(BandS)).c_str(),"[0]",0,1000);
1587     JZBratio[BandS]->Fit(pol0[BandS],"Q0R","",0,30);
1588    
1589    
1590    
1591    
1592    
1593    
1594     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("Bpred_mc"+any2string(canvasCounter)).c_str());
1595     can[canvasCounter]->cd();
1596     can[canvasCounter]->SetLogy(1);
1597     canvasCounter++;
1598     jzbprime_ossf_JZBP[allB]->SetMarkerStyle(4);
1599     jzbprime_ossf_JZBP[allB]->GetXaxis()->SetRangeUser(0,JZBmaxMC);
1600     jzbprime_ossf_JZBP[allB]->Draw("e1");
1601     jzbprime_ossf_JZBP[allB]->SetMinimum(0.0001);
1602     jzbprime_ossf_JZBP[allB]->SetMaximum(1000);
1603     Bpred[allB]->SetLineColor(kRed);
1604     Bpred[allB]->SetLineWidth(1);
1605     Bpred[allB]->Draw("hist,same");
1606     BpredSysP[allB]->SetLineColor(kBlack);
1607     BpredSysP[allB]->SetLineStyle(1);
1608     if(plotSys)BpredSysP[allB]->Draw("hist,same");
1609     BpredSysN[allB]->SetLineColor(kBlue);
1610     BpredSysN[allB]->SetLineStyle(1);
1611     if(plotSys)BpredSysN[allB]->Draw("hist,same");
1612     paveTitleMC->Draw("same");
1613     jzbprime_ossf_JZBP[allB]->Draw("e1,same");
1614     leg[canvasCounter] = new TLegend(0.60,0.64,0.94,0.77);
1615     leg[canvasCounter] ->SetFillStyle(4000);
1616     leg[canvasCounter] ->SetFillColor(kWhite);
1617     leg[canvasCounter] ->SetTextFont(42);
1618     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[allB],"MC true B","p");
1619     leg[canvasCounter] ->AddEntry(Bpred[allB],"data-driven B","l");
1620     leg[canvasCounter] ->Draw("same");
1621     eventSelectionPaveText->Draw("same");
1622     can[canvasCounter-1]->Update();
1623     saveCanvas(canvasCounter-1,"Bpred_mc");
1624    
1625    
1626     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_ratio_mc"+any2string(canvasCounter)).c_str());
1627     can[canvasCounter]->cd();
1628     canvasCounter++;
1629     float JZBzoom = 350;
1630     JZBratio[allB]->GetXaxis()->SetRangeUser(0,JZBzoom);
1631     JZBratio[allB]->GetYaxis()->SetTitle("Observed/Predicted");
1632     JZBratio[allB]->GetXaxis()->SetTitle("JZB (GeV)");
1633     JZBratio[allB]->GetYaxis()->SetNdivisions(519);
1634     JZBratio[allB]->GetYaxis()->SetRangeUser(-0.5,5.0);
1635     JZBratio[allB]->Draw("AP");
1636     paveTitleMC->Draw("same");
1637     paveText[canvasCounter-1] = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
1638     paveText[canvasCounter-1]->SetFillStyle(4000);
1639     paveText[canvasCounter-1]->SetFillColor(kWhite);
1640     paveText[canvasCounter-1]->SetTextFont(42);
1641     ostringstream jzb_agreement_allB_text;
1642     jzb_agreement_allB_text<< setprecision(2) << "mean =" << pol0[allB]->GetParameter(0) << " #pm " << setprecision(1) << pol0[allB]->GetParError(0);
1643     paveText[canvasCounter-1]->AddText("MC closure test");
1644     paveText[canvasCounter-1]->AddText(jzb_agreement_allB_text.str().c_str());
1645     paveText[canvasCounter-1]->Draw("same");
1646     pol0[allB]->Draw("same");
1647     f1line[canvasCounter][1] = new TF1("","1.2",0,1000);
1648     f1line[canvasCounter][2] = new TF1("","0.6",0,1000);
1649     f1line[canvasCounter][1]->SetLineColor(kBlue);
1650     f1line[canvasCounter][1]->SetLineStyle(2);
1651     f1line[canvasCounter][2]->SetLineColor(kBlue);
1652     f1line[canvasCounter][2]->SetLineStyle(2);
1653     f1line[canvasCounter][1]->Draw("same");
1654     f1line[canvasCounter][2]->Draw("same");
1655     f1line[canvasCounter][3] = new TF1("","1.0",0,1000);
1656     f1line[canvasCounter][3]->SetLineColor(kBlue);
1657     f1line[canvasCounter][3]->SetLineStyle(1);
1658     f1line[canvasCounter][3]->Draw("same");
1659     eventSelectionPaveText->Draw("same");
1660     eventSelectionPaveText->Draw("same");
1661     leg[canvasCounter] = new TLegend(0.60,0.5,0.94,0.77);
1662     leg[canvasCounter] ->SetFillStyle(4000);
1663     leg[canvasCounter] ->SetFillColor(kWhite);
1664     leg[canvasCounter] ->SetTextFont(42);
1665     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][1],"+20\% sys envelope","l");
1666     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][3],"ratio = 1","l");
1667     leg[canvasCounter] ->AddEntry(pol0[allB],"fit in [0,30] GeV","l");
1668     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][2],"-40\% sys envelope","l");
1669     leg[canvasCounter] ->Draw("same");
1670     can[canvasCounter-1]->Update();
1671     saveCanvas(canvasCounter-1,"jzb_ratio_mc");
1672    
1673     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_ratio_unzoomed_mc"+any2string(canvasCounter)).c_str());
1674     can[canvasCounter]->cd();
1675     canvasCounter++;
1676     JZBratio[allB]->GetXaxis()->SetRangeUser(0,600);
1677     JZBratio[allB]->GetYaxis()->SetTitle("Observed/Predicted");
1678     JZBratio[allB]->GetXaxis()->SetTitle("JZB (GeV)");
1679     JZBratio[allB]->GetYaxis()->SetNdivisions(519);
1680     JZBratio[allB]->GetYaxis()->SetRangeUser(0.0,9.0);
1681     JZBratio[allB]->Draw("AP");
1682     paveTitleMC->Draw("same");
1683     paveText[canvasCounter-1] = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
1684     paveText[canvasCounter-1]->SetFillStyle(4000);
1685     paveText[canvasCounter-1]->SetFillColor(kWhite);
1686     paveText[canvasCounter-1]->SetTextFont(42);
1687     ostringstream jzb_agreement_allB_unzoomed_text;
1688     jzb_agreement_allB_unzoomed_text<< setprecision(2) << "mean =" << pol0[allB]->GetParameter(0) << " #pm " << setprecision(1) << pol0[allB]->GetParError(0);
1689     paveText[canvasCounter-1]->AddText("MC closure test");
1690     paveText[canvasCounter-1]->AddText(jzb_agreement_allB_unzoomed_text.str().c_str());
1691     paveText[canvasCounter-1]->Draw("same");
1692     pol0[allB]->Draw("same");
1693     f1line[canvasCounter][1] = new TF1("","1.2",0,1000);
1694     f1line[canvasCounter][2] = new TF1("","0.6",0,1000);
1695     f1line[canvasCounter][1]->SetLineColor(kBlue);
1696     f1line[canvasCounter][1]->SetLineStyle(2);
1697     f1line[canvasCounter][2]->SetLineColor(kBlue);
1698     f1line[canvasCounter][2]->SetLineStyle(2);
1699     f1line[canvasCounter][1]->Draw("same");
1700     f1line[canvasCounter][2]->Draw("same");
1701     f1line[canvasCounter][3] = new TF1("","1.0",0,1000);
1702     f1line[canvasCounter][3]->SetLineColor(kBlue);
1703     f1line[canvasCounter][3]->SetLineStyle(1);
1704     f1line[canvasCounter][3]->Draw("same");
1705     eventSelectionPaveText->Draw("same");
1706     eventSelectionPaveText->Draw("same");
1707     leg[canvasCounter] = new TLegend(0.56,0.32,0.90,0.60);
1708     leg[canvasCounter] ->SetFillStyle(4000);
1709     leg[canvasCounter] ->SetFillColor(kWhite);
1710     leg[canvasCounter] ->SetTextFont(42);
1711     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][1],"+20\% sys envelope","l");
1712     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][3],"ratio = 1","l");
1713     leg[canvasCounter] ->AddEntry(pol0[allB],"fit in [0,30] GeV","l");
1714     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][2],"-40\% sys envelope","l");
1715     leg[canvasCounter] ->Draw("same");
1716     can[canvasCounter-1]->Update();
1717     saveCanvas(canvasCounter-1,"jzb_ratio_unzoomed_mc");
1718    
1719    
1720     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_ratio_BandS_mc"+any2string(canvasCounter)).c_str());
1721     can[canvasCounter]->cd();
1722     canvasCounter++;
1723     JZBratio[BandS]->GetXaxis()->SetRangeUser(0,600);
1724     JZBratio[BandS]->GetYaxis()->SetTitle("Observed/Predicted");
1725     JZBratio[BandS]->GetXaxis()->SetTitle("JZB (GeV)");
1726     JZBratio[BandS]->GetYaxis()->SetNdivisions(519);
1727     JZBratio[BandS]->GetYaxis()->SetRangeUser(0.0,9.0);
1728     JZBratio[BandS]->Draw("AP");
1729     paveTitleMC->Draw("same");
1730     paveText[canvasCounter-1] = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
1731     paveText[canvasCounter-1]->SetFillStyle(4000);
1732     paveText[canvasCounter-1]->SetFillColor(kWhite);
1733     paveText[canvasCounter-1]->SetTextFont(42);
1734     ostringstream jzb_agreement_BandS_text;
1735     jzb_agreement_BandS_text<< setprecision(2) << "mean =" << pol0[BandS]->GetParameter(0) << " #pm " << setprecision(1) << pol0[BandS]->GetParError(0);
1736     paveText[canvasCounter-1]->AddText("MC B+S closure test");
1737     paveText[canvasCounter-1]->AddText(jzb_agreement_BandS_text.str().c_str());
1738     paveText[canvasCounter-1]->Draw("same");
1739     pol0[BandS]->Draw("same");
1740     f1line[canvasCounter][1] = new TF1("","1.2",0,1000);
1741     f1line[canvasCounter][2] = new TF1("","0.6",0,1000);
1742     f1line[canvasCounter][1]->SetLineColor(kBlue);
1743     f1line[canvasCounter][1]->SetLineStyle(2);
1744     f1line[canvasCounter][2]->SetLineColor(kBlue);
1745     f1line[canvasCounter][2]->SetLineStyle(2);
1746     f1line[canvasCounter][1]->Draw("same");
1747     f1line[canvasCounter][2]->Draw("same");
1748     f1line[canvasCounter][3] = new TF1("","1.0",0,1000);
1749     f1line[canvasCounter][3]->SetLineColor(kBlue);
1750     f1line[canvasCounter][3]->SetLineStyle(1);
1751     f1line[canvasCounter][3]->Draw("same");
1752     eventSelectionPaveText->Draw("same");
1753     eventSelectionPaveText->Draw("same");
1754     leg[canvasCounter] = new TLegend(0.56,0.32,0.90,0.60);
1755     leg[canvasCounter] ->SetFillStyle(4000);
1756     leg[canvasCounter] ->SetFillColor(kWhite);
1757     leg[canvasCounter] ->SetTextFont(42);
1758     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][1],"+20\% sys envelope","l");
1759     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][3],"ratio = 1","l");
1760     leg[canvasCounter] ->AddEntry(pol0[allB],"fit in [0,30] GeV","l");
1761     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][2],"-40\% sys envelope","l");
1762     leg[canvasCounter] ->Draw("same");
1763     can[canvasCounter-1]->Update();
1764     saveCanvas(canvasCounter-1,"jzb_ratio_BandS_mc");
1765    
1766    
1767    
1768     float rData=1.00;
1769     Bpred[data] = (TH1F*)jzbprime_ossf_JZBN[data]->Clone();
1770     Bpred[data]->Add(jzbprime_osof_JZBN[data],-rData);
1771     Bpred[data]->Add(jzbprime_osof_JZBP[data],+rData);
1772    
1773     Bpredfin[data] = (TH1F*)jzbprime_ossf_JZBNfin[data]->Clone();
1774     Bpredfin[data]->Add(jzbprime_osof_JZBNfin[data],-rData);
1775     Bpredfin[data]->Add(jzbprime_osof_JZBPfin[data],+rData);
1776    
1777     float rDataSysP=1.00;
1778     BpredSysP[data] = (TH1F*)jzbprime_ossf_JZBN_MisCalibP[data]->Clone();
1779     BpredSysP[data]->Add(jzbprime_osof_JZBN_MisCalibP[data],-rDataSysP);
1780     BpredSysP[data]->Add(jzbprime_osof_JZBP_MisCalibP[data],+rDataSysP);
1781    
1782     float rDataSysN=1.00;
1783     BpredSysN[data] = (TH1F*)jzbprime_ossf_JZBN_MisCalibN[data]->Clone();
1784     BpredSysN[data]->Add(jzbprime_osof_JZBN_MisCalibN[data],-rDataSysN);
1785     BpredSysN[data]->Add(jzbprime_osof_JZBP_MisCalibN[data],+rDataSysN);
1786    
1787    
1788    
1789    
1790     // #--- fit Bpred[data]
1791     cb[data] = new TF1("cb12",InvCrystalBall,0,1000,5);
1792     cb[data]->SetParameter(0,Bpred[data]->GetBinContent(1));
1793     cb[data]->SetParameter(1,0.);
1794     cb[data]->SetParameter(2,f1Data->GetParameter(2));
1795     cb[data]->SetParameter(3,1.8);
1796     cb[data]->SetParameter(4,2.5);
1797     TH1F *Bpredrebindata=(TH1F*)Bpred[data]->Clone();
1798     Bpredrebindata->Rebin();
1799     Bpredrebindata->Fit(cb[data],"N0");
1800     // Bpred[data]->Fit(cb[data],"N0");
1801     cb[data]->SetLineColor(kBlue);
1802     cout << "chi2 = "<<cb[data]->GetChisquare()<<" ndf = "<<cb[data]->GetNDF()<< " chi2/ndf = "<< cb[data]->GetChisquare()/cb[data]->GetNDF()<<endl;
1803    
1804     cbP[data] = new TF1("cbP12",InvCrystalBallP,0,1000,5);
1805     cbP[data]->SetParameter(0,cb[data]->GetParameter(0));
1806     cbP[data]->SetParameter(1,cb[data]->GetParameter(1));
1807     cbP[data]->SetParameter(2,cb[data]->GetParameter(2));
1808     cbP[data]->SetParameter(3,cb[data]->GetParameter(3));
1809     cbP[data]->SetParameter(4,cb[data]->GetParameter(4));
1810     cbP[data]->SetLineStyle(2);
1811     cbP[data]->SetLineColor(kBlue);
1812    
1813    
1814     cbN[data] = new TF1("cbN12",InvCrystalBallN,0,1000,5);
1815     cbN[data]->SetParameter(0,cb[data]->GetParameter(0));
1816     cbN[data]->SetParameter(1,cb[data]->GetParameter(1));
1817     cbN[data]->SetParameter(2,cb[data]->GetParameter(2));
1818     cbN[data]->SetParameter(3,cb[data]->GetParameter(3));
1819     cbN[data]->SetParameter(4,cb[data]->GetParameter(4));
1820     cbN[data]->SetLineStyle(2);
1821     cbN[data]->SetLineColor(kBlue);
1822    
1823     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("Bpred_data"+any2string(canvasCounter)).c_str());
1824     can[canvasCounter]->cd();
1825     can[canvasCounter]->SetLogy(1);
1826     canvasCounter++;
1827     TH1F *Bpred_data_pres=(TH1F*)Bpred[data]->Clone();
1828     Bpred_data_pres->Rebin();
1829     TH1F *jzbprime_ossfJZBPpresent=(TH1F*)jzbprime_ossf_JZBP[data]->Clone();
1830     jzbprime_ossfJZBPpresent->Rebin();
1831     // jzbprime_ossf_JZBP[data]->GetXaxis()->SetRangeUser(0,350);
1832     jzbprime_ossfJZBPpresent->GetXaxis()->SetRangeUser(0,350);
1833     // jzbprime_ossf_JZBP[data]->Draw("e1");
1834     jzbprime_ossfJZBPpresent->Draw("e1");
1835     // Bpred[data]->SetLineColor(kRed);
1836     Bpred_data_pres->SetLineColor(kRed);
1837     // Bpred[data]->SetLineStyle(0);
1838     Bpred_data_pres->SetLineStyle(0);
1839     // Bpred[data]->Draw("hist,same");
1840     Bpred_data_pres->Draw("hist,same");
1841     // cb[data]->Draw("same");cbN[data]->Draw("same");cbP[data]->Draw("same");
1842     paveTitle->Draw("same");
1843     jzbprime_ossfJZBPpresent->Draw("e1,same");
1844     jzbprime_ossfJZBPpresent->Draw("axis,same");
1845     // jzbprime_ossf_JZBP[data]->Draw("e1,same");
1846     // jzbprime_ossf_JZBP[data]->Draw("axis,same");
1847     cb[data]->Draw("same");cbN[data]->Draw("same");cbP[data]->Draw("same");
1848     leg[canvasCounter] = new TLegend(0.60,0.5,0.94,0.77);
1849     leg[canvasCounter] ->SetFillStyle(4000);
1850     leg[canvasCounter] ->SetFillColor(kWhite);
1851     leg[canvasCounter] ->SetTextFont(42);
1852     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[data],"observed (data)","p");
1853     leg[canvasCounter] ->AddEntry(Bpred[data],"predicted (data)","l");
1854     leg[canvasCounter] ->AddEntry(cb[data],"predicted fit (data)","l");
1855     leg[canvasCounter] ->AddEntry(cbP[data],"stat. uncert.","l");
1856     leg[canvasCounter] ->Draw("same");
1857     eventSelectionPaveText->Draw("same");
1858     can[canvasCounter-1]->Update();
1859     saveCanvas(canvasCounter-1,"Bpred_data");
1860     /*
1861     cout << "hello" << endl; // MAB HACK
1862     cout << "JZB Prime OSSF JZBP (data) : " << endl;
1863     for(int i=1;i<=jzbprime_ossf_JZBP[data]->GetNbinsX();i++)
1864     {
1865     cout << i << "\t" << jzbprime_ossf_JZBP[data]->GetBinCenter(i) << "\t" << jzbprime_ossf_JZBP[data]->GetBinContent(i) << endl;
1866     }
1867     cout << "Bpred[data] : " << endl;
1868     for(int i=1;i<=Bpred[data]->GetNbinsX();i++)
1869     {
1870     cout << i << "\t" << Bpred[data]->GetBinCenter(i) << "\t" << Bpred[data]->GetBinContent(i) << endl;
1871     }
1872     */
1873    
1874    
1875    
1876    
1877    
1878    
1879     JZBratio[data] = histRatio(jzbprime_ossf_JZBP[data], Bpred[data], data);
1880     pol0[data] = new TF1(("s"+any2string(data)).c_str(),"[0]",0,1000);
1881     JZBratio[data]->Fit(pol0[data],"Q0R","",0,30);
1882     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_ratio_data"+any2string(canvasCounter)).c_str());
1883     can[canvasCounter]->cd();
1884     canvasCounter++;
1885     JZBratio[data]->GetXaxis()->SetRangeUser(0,350);
1886     JZBratio[data]->GetYaxis()->SetTitle("Observed/Predicted");
1887     JZBratio[data]->GetXaxis()->SetTitle("JZB (GeV)");
1888     JZBratio[data]->GetYaxis()->SetNdivisions(519);
1889     JZBratio[data]->GetYaxis()->SetRangeUser(-0.5,5.0);
1890     JZBratio[data]->Draw("AP");
1891     paveTitle->Draw("same");
1892     paveText[canvasCounter-1] = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
1893     paveText[canvasCounter-1]->SetFillStyle(4000);
1894     paveText[canvasCounter-1]->SetFillColor(kWhite);
1895     paveText[canvasCounter-1]->SetTextFont(42);
1896     ostringstream jzb_agreement_data_text;
1897     jzb_agreement_data_text<< setprecision(2) << "mean =" << pol0[data]->GetParameter(0) << " #pm " << setprecision(1) << pol0[data]->GetParError(0);
1898     paveText[canvasCounter-1]->AddText("Data closure test");
1899     paveText[canvasCounter-1]->AddText(jzb_agreement_data_text.str().c_str());
1900     paveText[canvasCounter-1]->Draw("same");
1901     pol0[data]->Draw("same");
1902     f1line[canvasCounter][1] = new TF1("","1.2",0,1000);
1903     f1line[canvasCounter][2] = new TF1("","0.6",0,1000);
1904     f1line[canvasCounter][1]->SetLineColor(kBlue);
1905     f1line[canvasCounter][1]->SetLineStyle(2);
1906     f1line[canvasCounter][2]->SetLineColor(kBlue);
1907     f1line[canvasCounter][2]->SetLineStyle(2);
1908     f1line[canvasCounter][1]->Draw("same");
1909     f1line[canvasCounter][2]->Draw("same");
1910     f1line[canvasCounter][3] = new TF1("","1.0",0,1000);
1911     f1line[canvasCounter][3]->SetLineColor(kBlue);
1912     f1line[canvasCounter][3]->SetLineStyle(1);
1913     f1line[canvasCounter][3]->Draw("same");
1914     leg[canvasCounter] = new TLegend(0.60,0.5,0.94,0.77);
1915     leg[canvasCounter] ->SetFillStyle(4000);
1916     leg[canvasCounter] ->SetFillColor(kWhite);
1917     leg[canvasCounter] ->SetTextFont(42);
1918     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][1],"+20\% sys envelope","l");
1919     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][3],"ratio = 1","l");
1920     leg[canvasCounter] ->AddEntry(pol0[data],"fit in [0,30] GeV","l");
1921     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][2],"-40\% sys envelope","l");
1922     leg[canvasCounter] ->Draw("same");
1923     // eventSelectionPaveText->Draw("same");
1924     eventSelectionPaveText->Draw("same");
1925     can[canvasCounter-1]->Update();
1926     // saveCanvas(canvasCounter-1,"jzb_ratio_data");
1927    
1928     //********* rebinned one below!!!!!
1929    
1930    
1931    
1932     // JZBratiofin[data] = (TH1F*)jzbprime_ossf_JZBPfin[data]->Clone();
1933     // JZBratiofin[data]->Divide(Bpredfin[data]);
1934     JZBratiofin[data] = histRatiofin(jzbprime_ossf_JZBPfin[data], Bpredfin[data], data);
1935     // histRatio(jzbprime_ossf_JZBPfin[data], Bpredfin[data], data);
1936     pol0[12] = new TF1(("s"+any2string(data)).c_str(),"[0]",0,1000);
1937     JZBratiofin[data]->Fit(pol0[12],"Q0R","",0,30);
1938     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_ratio_data"+any2string(canvasCounter)).c_str());
1939     can[canvasCounter]->cd();
1940     canvasCounter++;
1941     JZBratiofin[data]->GetXaxis()->SetRangeUser(0,350);
1942     JZBratiofin[data]->GetYaxis()->SetTitle("Observed/Predicted");
1943     JZBratiofin[data]->GetXaxis()->SetTitle("JZB (GeV)");
1944     JZBratiofin[data]->GetYaxis()->SetNdivisions(519);
1945     JZBratiofin[data]->GetYaxis()->SetRangeUser(-0.5,5.0);
1946     JZBratiofin[data]->Draw("AP");
1947     paveTitle->Draw("same");
1948     paveText[canvasCounter-1] = new TPaveText(0.2,0.78,0.49,0.91,"blNDC");
1949     paveText[canvasCounter-1]->SetFillStyle(4000);
1950     paveText[canvasCounter-1]->SetFillColor(kWhite);
1951     paveText[canvasCounter-1]->SetTextFont(42);
1952     ostringstream jzbf_agreement_data_text;
1953     jzbf_agreement_data_text<< setprecision(2) << "mean =" << pol0[12]->GetParameter(0) << " #pm " << setprecision(1) << pol0[12]->GetParError(0);
1954     paveText[canvasCounter-1]->AddText("Data closure test");
1955     paveText[canvasCounter-1]->AddText(jzb_agreement_data_text.str().c_str());
1956     paveText[canvasCounter-1]->Draw("same");
1957     pol0[12]->Draw("same");
1958     f1line[canvasCounter][1] = new TF1("","1.2",0,1000);
1959     f1line[canvasCounter][2] = new TF1("","0.6",0,1000);
1960     f1line[canvasCounter][1]->SetLineColor(kBlue);
1961     f1line[canvasCounter][1]->SetLineStyle(2);
1962     f1line[canvasCounter][2]->SetLineColor(kBlue);
1963     f1line[canvasCounter][2]->SetLineStyle(2);
1964     f1line[canvasCounter][1]->Draw("same");
1965     f1line[canvasCounter][2]->Draw("same");
1966     f1line[canvasCounter][3] = new TF1("","1.0",0,1000);
1967     f1line[canvasCounter][3]->SetLineColor(kBlue);
1968     f1line[canvasCounter][3]->SetLineStyle(1);
1969     f1line[canvasCounter][3]->Draw("same");
1970     leg[canvasCounter] = new TLegend(0.60,0.5,0.94,0.77);
1971     leg[canvasCounter] ->SetFillStyle(4000);
1972     leg[canvasCounter] ->SetFillColor(kWhite);
1973     leg[canvasCounter] ->SetTextFont(42);
1974     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][1],"+20\% sys envelope","l");
1975     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][3],"ratio = 1","l");
1976     leg[canvasCounter] ->AddEntry(pol0[12],"fit in [0,30] GeV","l");
1977     leg[canvasCounter] ->AddEntry(f1line[canvasCounter][2],"-40\% sys envelope","l");
1978     leg[canvasCounter] ->Draw("same");
1979     // eventSelectionPaveText->Draw("same");
1980     eventSelectionPaveText->Draw("same");
1981     can[canvasCounter-1]->Update();
1982     saveCanvas(canvasCounter-1,"jzb_ratio_data_REBINNED");
1983    
1984    
1985    
1986    
1987     //********* rebinned one above!
1988     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_zjet_vs_ttbar_vs_lm4_mc"+any2string(canvasCounter)).c_str());
1989     can[canvasCounter]->cd();
1990     can[canvasCounter]->SetLogy(1);
1991     canvasCounter++;
1992     jzb_ossf_rebinned[zjets] = (TH1F*)jzb_ossf[zjets]->Clone();
1993     jzb_ossf_rebinned[ttbar] = (TH1F*)jzb_ossf[ttbar]->Clone();
1994     jzb_ossf_rebinned[lm4] = (TH1F*)jzb_ossf[lm4]->Clone();
1995     jzb_ossf_rebinned[zjets]->Rebin(4);
1996     jzb_ossf_rebinned[ttbar]->Rebin(4);
1997     jzb_ossf_rebinned[lm4]->Rebin(4);
1998     jzb_ossf_rebinned[ttbar]->SetFillColor(0);
1999     jzb_ossf_rebinned[lm4]->SetFillColor(0);
2000     jzb_ossf_rebinned[ttbar]->SetLineColor(kMagenta+2);
2001     jzb_ossf_rebinned[lm4]->SetLineColor(kViolet+7);
2002     jzb_ossf_rebinned[lm4]->SetLineStyle(kViolet+7);
2003     jzb_ossf_rebinned[ttbar]->SetLineWidth(2);
2004     jzb_ossf_rebinned[lm4]->SetLineWidth(2);
2005     jzb_ossf_rebinned[lm4]->SetLineStyle(2);
2006     jzb_ossf_rebinned[zjets]->GetXaxis()->SetRangeUser(-200,400);
2007     jzb_ossf_rebinned[zjets]->SetMaximum(1000);
2008     jzb_ossf_rebinned[zjets]->SetMinimum(0.04);
2009     jzb_ossf_rebinned[zjets]->Draw("hist");
2010     jzb_ossf_rebinned[ttbar]->Draw("hist,same");
2011     jzb_ossf_rebinned[lm4]->Draw("hist,same");
2012     paveTitleMC->Draw("same");
2013     leg[canvasCounter] = new TLegend(0.60,0.64,0.94,0.77);
2014     leg[canvasCounter] ->SetFillStyle(4000);
2015     leg[canvasCounter] ->SetFillColor(kWhite);
2016     leg[canvasCounter] ->SetTextFont(42);
2017     leg[canvasCounter] ->AddEntry(jzb_ossf_rebinned[zjets],"ZJets","f");
2018     leg[canvasCounter] ->AddEntry(jzb_ossf_rebinned[ttbar],"ttbar","f");
2019     leg[canvasCounter] ->AddEntry(jzb_ossf_rebinned[lm4],"signal x 20","f");
2020     leg[canvasCounter] ->Draw("same");
2021     eventSelectionPaveText->Draw("same");
2022     can[canvasCounter-1]->Update();
2023     saveCanvas(canvasCounter-1,"jzb_zjet_vs_ttbar_vs_lm4_mc");
2024    
2025     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("ttbar_emu_mc"+any2string(canvasCounter)).c_str());
2026     can[canvasCounter]->cd();
2027     can[canvasCounter]->SetLogy(1);
2028     canvasCounter++;
2029     jzb_osof_rebinned[ttbar] = (TH1F*)jzb_osof[ttbar]->Clone();
2030     jzb_osof_rebinned[ttbar]->Rebin(4);
2031     jzb_osof_rebinned[ttbar]->SetFillColor(0);
2032     jzb_osof_rebinned[ttbar]->SetLineColor(kMagenta+2);
2033     jzb_osof_rebinned[ttbar]->SetLineWidth(2);
2034     jzb_osof_rebinned[ttbar]->SetLineStyle(2);
2035     jzb_osof_rebinned[ttbar]->GetXaxis()->SetRangeUser(-100,400);
2036     jzb_osof_rebinned[ttbar]->SetMinimum(0.001);
2037     jzb_osof_rebinned[ttbar]->SetMaximum(10);
2038     jzb_osof_rebinned[ttbar]->Draw("hist");
2039     jzb_ossf_rebinned[ttbar]->Draw("hist,same");
2040     paveTitleMC->Draw("same");
2041     leg[canvasCounter] = new TLegend(0.60,0.64,0.94,0.77);
2042     leg[canvasCounter] ->SetFillStyle(4000);
2043     leg[canvasCounter] ->SetFillColor(kWhite);
2044     leg[canvasCounter] ->SetTextFont(42);
2045     leg[canvasCounter] ->AddEntry(jzb_ossf_rebinned[ttbar],"ttbar #rightarrow (ee or #mu#mu)","f");
2046     leg[canvasCounter] ->AddEntry(jzb_osof_rebinned[ttbar],"ttbar #rightarrow e#mu ","f");
2047     leg[canvasCounter] ->Draw("same");
2048     eventSelectionPaveText->Draw("same");
2049     can[canvasCounter-1]->Update();
2050     saveCanvas(canvasCounter-1,"ttbar_emu_mc");
2051    
2052     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("mll_effratio_mc"+any2string(canvasCounter)).c_str());
2053     can[canvasCounter]->cd();
2054     can[canvasCounter]->SetLogy(1);
2055     canvasCounter++;
2056     mll_mumu[zjets]->SetFillColor(0);
2057     mll_mumu[zjets]->SetLineColor(kBlue);
2058     mll_mumu[zjets]->SetLineWidth(2);
2059     mll_mumu[zjets]->SetLineStyle(1);
2060     mll_ee[zjets]->SetFillColor(0);
2061     mll_ee[zjets]->SetLineColor(kOrange+2);
2062     mll_ee[zjets]->SetLineWidth(2);
2063     mll_ee[zjets]->SetLineStyle(2);
2064     mll_mumu[zjets]->GetXaxis()->SetRangeUser(50,150);
2065     mll_mumu[zjets]->SetMaximum(1500);
2066     mll_mumu[zjets]->SetMinimum(0.1);
2067     mll_mumu[zjets]->Draw("hist");
2068     mll_ee[zjets]->Draw("hist,same");
2069     paveTitleMC->Draw("same");
2070     leg[canvasCounter] = new TLegend(0.60,0.64,0.94,0.77);
2071     leg[canvasCounter] ->SetFillStyle(4000);
2072     leg[canvasCounter] ->SetFillColor(kWhite);
2073     leg[canvasCounter] ->SetTextFont(42);
2074     leg[canvasCounter] ->AddEntry(mll_mumu[zjets],"Zjets #rightarrow #mu#mu","f");
2075     leg[canvasCounter] ->AddEntry(mll_ee[zjets],"Zjets #rightarrow ee ","f");
2076     leg[canvasCounter] ->Draw("same");
2077     line[canvasCounter][1] = new TLine(70. , .0,70. ,1.);
2078     line[canvasCounter][2] = new TLine(110., .0,110.,1.0);
2079     line[canvasCounter][1]->SetLineColor(kBlue);
2080     line[canvasCounter][1]->SetLineStyle(2);
2081     line[canvasCounter][2]->SetLineColor(kBlue);
2082     line[canvasCounter][2]->SetLineStyle(2);
2083     line[canvasCounter][1]->Draw("same");
2084     line[canvasCounter][2]->Draw("same");
2085     eventSelectionPaveText->Draw("same");
2086     eventSelectionPaveText->Draw("same");
2087     can[canvasCounter-1]->Update();
2088     saveCanvas(canvasCounter-1,"mll_effratio_mc");
2089    
2090     jzb_ossf_rebinned[data] = (TH1F*)jzb_ossf[data]->Clone() ; jzb_osof_rebinned[data] = (TH1F*)jzb_osof[data]->Clone();
2091     jzb_ossf_rebinned[data]->Rebin(4);
2092     jzb_osof_rebinned[data]->Rebin(4);
2093     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("jzb_eemumu_emu_data"+any2string(canvasCounter)).c_str());
2094     can[canvasCounter]->cd();
2095     can[canvasCounter]->SetLogy(1);
2096     canvasCounter++;
2097     jzb_ossf_rebinned[data]->SetMarkerStyle(4);
2098     jzb_osof_rebinned[data]->SetMarkerStyle(25);
2099     jzb_osof_rebinned[data]->SetLineColor(kRed);
2100     jzb_osof_rebinned[data]->SetMarkerColor(kRed);
2101     // jzb_ossf_rebinned[data]->GetXaxis()->SetRangeUser(-150,150);
2102     jzb_ossf_rebinned[data]->GetXaxis()->SetRangeUser(-120,350);
2103     jzb_ossf_rebinned[data]->Draw("e1");
2104     jzb_osof_rebinned[data]->Draw("e1,same");
2105     paveTitle->Draw("same");
2106     leg[canvasCounter] = new TLegend(0.68,0.60,0.93,0.78);
2107     leg[canvasCounter]->SetHeader("data");
2108     leg[canvasCounter] ->SetFillStyle(4000);
2109     leg[canvasCounter] ->SetFillColor(kWhite);
2110     leg[canvasCounter] ->SetTextFont(42);
2111     leg[canvasCounter] ->AddEntry(jzb_ossf_rebinned[data],"#mu#mu or ee","p");
2112     leg[canvasCounter] ->AddEntry(jzb_osof_rebinned[data],"e#mu","p");
2113     leg[canvasCounter] ->Draw("same");
2114     eventSelectionPaveText->Draw("same");
2115     can[canvasCounter-1]->Update();
2116     saveCanvas(canvasCounter-1,"jzb_eemumu_emu_data");
2117    
2118    
2119     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("BpredBandS_mc"+any2string(canvasCounter)).c_str());
2120     can[canvasCounter]->cd();
2121     can[canvasCounter]->SetLogy(1);
2122     canvasCounter++;
2123     jzbprime_ossf_JZBP[BandS]->SetMarkerStyle(4);
2124     jzbprime_ossf_JZBP[BandS]->GetXaxis()->SetRangeUser(0,300);
2125     jzbprime_ossf_JZBP[BandS]->Draw("e1");
2126     jzbprime_ossf_JZBP[BandS]->SetMinimum(0.005);
2127     jzbprime_ossf_JZBP[BandS]->SetMaximum(1000);
2128     jzbprime_ossf_JZBP[allB]->SetLineStyle(3);
2129     jzbprime_ossf_JZBP[allB]->Draw("hist,same");
2130     Bpred[BandS]->SetLineColor(kRed);
2131     Bpred[BandS]->SetLineWidth(1);
2132     Bpred[BandS]->Draw("hist,same");
2133     // BpredSysP[BandS]->SetLineColor(kBlack);
2134     // BpredSysP[BandS]->SetLineStyle(1);
2135     // if(plotSys)BpredSysP[BandS]->Draw("hist,same");
2136     // BpredSysN[BandS]->SetLineColor(kBlue);
2137     // BpredSysN[BandS]->SetLineStyle(1);
2138     // if(plotSys)BpredSysN[BandS]->Draw("hist,same");
2139     paveTitleMC->Draw("same");
2140     jzbprime_ossf_JZBP[BandS]->Draw("e1,same");
2141     jzbprime_ossf_JZBP[lm4]->Draw("hist,same");
2142     jzbprime_ossf_JZBP[lm4]->SetLineWidth(2);
2143     jzbprime_ossf_JZBP[lm4]->SetLineStyle(2);
2144     jzbprime_ossf_JZBP[lm4]->SetLineColor(kViolet+7);
2145     leg[canvasCounter] = new TLegend(0.63,0.57,0.91,0.77);
2146     leg[canvasCounter] ->SetFillStyle(4000);
2147     leg[canvasCounter] ->SetFillColor(kWhite);
2148     leg[canvasCounter] ->SetTextFont(42);
2149     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[BandS],"MC B + S","p");
2150     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[lm4],"signal #times 20","f");
2151     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[allB],"MC B","f");
2152     leg[canvasCounter] ->AddEntry(Bpred[BandS],"data-driven B","f");
2153     leg[canvasCounter] ->Draw("same");
2154     jzbprime_ossf_JZBP[BandS]->Draw("e1,same");
2155     jzbprime_ossf_JZBP[BandS]->Draw("axis,same");
2156     eventSelectionPaveText->Draw("same");
2157     can[canvasCounter-1]->Update();
2158     saveCanvas(canvasCounter-1,"BpredBandS_mc");
2159    
2160    
2161    
2162    
2163     // #--- fit Bpred[data]
2164     cb[zjets_un] = new TF1("cb0",InvCrystalBall,0,1000,5);
2165     cb[zjets_un]->SetParameter(0,jzbprime_ossf_JZBN[zjets_un]->GetBinContent(1));
2166     cb[zjets_un]->SetParameter(1,0.);
2167     cb[zjets_un]->SetParameter(2,f1MC->GetParameter(2));
2168     cb[zjets_un]->SetParameter(3,1.8);
2169     cb[zjets_un]->SetParameter(4,2.5);
2170     jzbprime_ossf_JZBN[zjets_un]->Fit(cb[zjets_un],"N0");
2171     cb[zjets_un]->SetLineColor(kBlue);
2172     cout << "chi2 = "<<cb[zjets_un]->GetChisquare()<<" ndf = "<<cb[zjets_un]->GetNDF()<< " chi2/ndf = "<< cb[zjets_un]->GetChisquare()/cb[zjets_un]->GetNDF()<<endl;
2173    
2174     cbP[zjets_un] = new TF1("cbP0",InvCrystalBallP,0,1000,5);
2175     cbP[zjets_un]->SetParameter(0,cb[zjets_un]->GetParameter(0));
2176     cbP[zjets_un]->SetParameter(1,cb[zjets_un]->GetParameter(1));
2177     cbP[zjets_un]->SetParameter(2,cb[zjets_un]->GetParameter(2));
2178     cbP[zjets_un]->SetParameter(3,cb[zjets_un]->GetParameter(3));
2179     cbP[zjets_un]->SetParameter(4,cb[zjets_un]->GetParameter(4));
2180     cbP[zjets_un]->SetLineStyle(2);
2181     cbP[zjets_un]->SetLineColor(kBlue);
2182    
2183    
2184     cbN[zjets_un] = new TF1("cbN0",InvCrystalBallN,0,1000,5);
2185     cbN[zjets_un]->SetParameter(0,cb[zjets_un]->GetParameter(0));
2186     cbN[zjets_un]->SetParameter(1,cb[zjets_un]->GetParameter(1));
2187     cbN[zjets_un]->SetParameter(2,cb[zjets_un]->GetParameter(2));
2188     cbN[zjets_un]->SetParameter(3,cb[zjets_un]->GetParameter(3));
2189     cbN[zjets_un]->SetParameter(4,cb[zjets_un]->GetParameter(4));
2190     cbN[zjets_un]->SetLineStyle(2);
2191     cbN[zjets_un]->SetLineColor(kBlue);
2192     can[canvasCounter] = new TCanvas(("c"+any2string(canvasCounter)).c_str(),("JZBN_vs_JZBP_ZJets_mc"+any2string(canvasCounter)).c_str());
2193     can[canvasCounter]->cd();
2194     can[canvasCounter]->SetLogy(1);
2195     canvasCounter++;
2196     jzbprime_ossf_JZBP[zjets_un]->GetXaxis()->SetRangeUser(0,120);
2197     jzbprime_ossf_JZBP[zjets_un]->SetMarkerStyle(24);
2198     jzbprime_ossf_JZBN[zjets_un]->SetLineColor(kRed);
2199     jzbprime_ossf_JZBP[zjets_un]->SetMaximum(5000);
2200     jzbprime_ossf_JZBP[zjets_un]->SetMinimum(0.1);
2201     jzbprime_ossf_JZBP[zjets_un]->Draw("e1");
2202     jzbprime_ossf_JZBN[zjets_un]->Draw("hist,same");
2203     cb[zjets_un]->Draw("same");
2204     cbP[zjets_un]->Draw("same");
2205     cbN[zjets_un]->Draw("same");
2206     paveTitleMCZjets->Draw("same");
2207     leg[canvasCounter] = new TLegend(0.60,0.5,0.94,0.77);
2208     leg[canvasCounter] ->SetFillStyle(4000);
2209     leg[canvasCounter] ->SetFillColor(kWhite);
2210     leg[canvasCounter] ->SetTextFont(42);
2211     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBP[zjets_un],"JZB>0","p");
2212     leg[canvasCounter] ->AddEntry(jzbprime_ossf_JZBN[zjets_un],"JZB<0","l");
2213     leg[canvasCounter] ->AddEntry(cb[zjets_un],"JZB<0 fit","l");
2214     leg[canvasCounter] ->AddEntry(cbP[zjets_un],"stat. uncert.","l");
2215     leg[canvasCounter] ->Draw("same");
2216     eventSelectionPaveText->Draw("same");
2217     can[canvasCounter-1]->Update();
2218     saveCanvas(canvasCounter-1,"JZBN_vs_JZBP_ZJets_mc");
2219    
2220     TCanvas *wingfinder = new TCanvas("wingfinder","Wing Finder");
2221     wingfinder->SetLogy(1);
2222     jzb_ossfee[data]->SetLineColor(kBlue);
2223     jzb_ossfee[data]->GetXaxis()->SetRangeUser(-100,100);
2224     jzb_ossfmm[data]->SetLineColor(kRed);
2225     jzb_ossfee[data]->GetYaxis()->SetTitle("Prob");
2226     jzb_ossfmm[data]->GetYaxis()->SetTitle("Prob");
2227     jzb_ossfee[data]->DrawNormalized("hist");
2228     jzb_ossfmm[data]->DrawNormalized("same,hist");
2229     leg[canvasCounter] = new TLegend(0.60,0.6,0.94,0.77);
2230     leg[canvasCounter] ->SetFillStyle(4000);
2231     leg[canvasCounter] ->SetFillColor(kWhite);
2232     leg[canvasCounter] ->SetTextFont(42);
2233     leg[canvasCounter] ->AddEntry(jzb_ossfee[data],"ee","l");
2234     leg[canvasCounter] ->AddEntry(jzb_ossfmm[data],"#mu#mu","l");
2235     leg[canvasCounter] ->Draw("same");
2236     eventSelectionPaveText->Draw("same");
2237     paveTitle->Draw("same");
2238     can[canvasCounter-1]->Update();
2239     saveCanvas(canvasCounter-1,"comparison_ee_mm");
2240    
2241    
2242     range = 50;
2243     vector<float> predvec;
2244     vector<float> obsvec;
2245     for(int i =0 ; i<6 ; i++)
2246     {
2247    
2248     float fraction = histIntegral(Bpred[data],range)/histIntegral(Bpred[data],0);
2249     float zjetsEstimate = histIntegral(jzbprime_ossf_JZBN[data],range);
2250     float zjetsEstimateMC = histIntegral(jzbprime_ossf_JZBN[allB],range);
2251     float ttbarEstimate = histIntegral(jzbprime_osof_JZBP[data],range) - histIntegral(jzbprime_osof_JZBN[data],range);
2252     float ttbarEstimateN = histIntegral(jzbprime_osof_JZBN[data],range);
2253     float ttbarEstimateP = histIntegral(jzbprime_osof_JZBP[data],range);
2254     // cout << "> "<< range << " : "<< histIntegral(jzbprime_ossf_JZBN[data],range) << " " << histIntegral(jzbprime_osof_JZBP[data],range) << " - " << histIntegral(jzbprime_osof_JZBN[data],range) <<endl;
2255     // cout<< histIntegral(jzbprime_ossf_JZBP[data],range) << " " << histIntegral(Bpred[data],range) << " >" << range << " : " << histIntegral(Bpred[data],range)/histIntegral(Bpred[data],0) <<endl;
2256    
2257     float obs = histIntegral(jzbprime_ossf_JZBP[data],range);
2258     float pred = histIntegral(jzbprime_ossf_JZBN[data],range) + histIntegral(jzbprime_osof_JZBP[data],range) - histIntegral(jzbprime_osof_JZBN[data],range);
2259     float predStatP = statErrorP(pred);
2260     float predStatN = statErrorN(pred);
2261    
2262    
2263    
2264     float mobs = histIntegral(jzbprime_ossf_JZBP[allB],range);
2265     float mpred = histIntegral(jzbprime_ossf_JZBN[allB],range) + histIntegral(jzbprime_osof_JZBP[allB],range) - histIntegral(jzbprime_osof_JZBN[allB],range);
2266     float mpredStatP = statErrorP(pred);
2267     float mpredStatN = statErrorN(pred);
2268    
2269     float allBSampleError = 0;
2270     float allBSample = histIntegralAndError(jzbprime_ossf_JZBP[allB],range,allBSampleError);
2271    
2272     /*
2273     float zjetsSampleError = 0;
2274     float zjetsSample = histIntegralAndError(jzbprime_ossf_JZBP[zjets],range,zjetsSampleError);
2275     float ttbarSampleError = 0;
2276     float ttbarSample = histIntegralAndError(jzbprime_ossf_JZBP[ttbar],range,ttbarSampleError);
2277     float lm4SampleError = 0;
2278     float lm4Sample = histIntegralAndError(jzbprime_ossf_JZBP[lm4],range,lm4SampleError);
2279     */
2280    
2281     // bin by bin variation
2282     float predSysP = histIntegral(jzbprime_ossf_JZBN_MisCalibP[data],range) + histIntegral(jzbprime_osof_JZBP_MisCalibP[data],range) - histIntegral(jzbprime_osof_JZBN_MisCalibP[data],range);
2283     float predSysN = histIntegral(jzbprime_ossf_JZBN_MisCalibN[data],range) + histIntegral(jzbprime_osof_JZBP_MisCalibN[data],range) - histIntegral(jzbprime_osof_JZBN_MisCalibN[data],range);
2284     float sysP = TMath::Abs(predSysP-pred);
2285     float sysN = TMath::Abs(pred - predSysN);
2286     sysN = sysN + pred*(1-pol0[allB]->GetParameter(0)); // consider sys bias
2287     //blublublu
2288    
2289     float mpredSysP = histIntegral(jzbprime_ossf_JZBN_MisCalibP[allB],range) + histIntegral(jzbprime_osof_JZBP_MisCalibP[allB],range) - histIntegral(jzbprime_osof_JZBN_MisCalibP[allB],range);
2290     float mpredSysN = histIntegral(jzbprime_ossf_JZBN_MisCalibN[allB],range) + histIntegral(jzbprime_osof_JZBP_MisCalibN[allB],range) - histIntegral(jzbprime_osof_JZBN_MisCalibN[allB],range);
2291     float msysP = TMath::Abs(mpredSysP-mpred);
2292     float msysN = TMath::Abs(mpred - mpredSysN);
2293     msysN = msysN + pred*(1-pol0[allB]->GetParameter(0)); // consider sys bias
2294    
2295    
2296     // cout << "> "<< range << " obs = " << obs << " pred = " << pred <<" (-" << pred*(1-pol0[allB]->GetParameter(0)) << " sys bias)";
2297     cout << "JZB > "<< range << " obs = " << obs << " pred = " << pred ;
2298     cout << " (sys) +" << setprecision(2) << sysP << "-" << sysN << " +- (stat) " << predStatP << "-" << predStatN << endl;
2299     cout << " ### zjetsEstimate = "<< zjetsEstimate << " ttbarEstimate = " << ttbarEstimate << " ttbarEstimateP = " << ttbarEstimateP << " ttbarEstimateN = " << ttbarEstimateN;
2300     cout << " f = " << setprecision(2) << fraction << "(%)" <<endl;
2301     cout << "MC allB = " << allBSample << " +- " << allBSampleError <<endl;
2302    
2303     cout << "MC: JZB > "<< range << " obs = " << mobs << " pred = " << mpred ;
2304     cout << " (sys) +" << setprecision(2) << msysP << "-" << msysN << " +- (stat) " << mpredStatP << "-" << mpredStatN << endl;
2305    
2306     predvec.push_back(pred);
2307     obsvec.push_back(obs);
2308     range += 10;
2309     }
2310     for(int i=0;i<predvec.size();i++) {
2311     cout << 10+i*10 << ";" << obsvec[i] << ";" << predvec[i] << endl;
2312     }
2313    
2314     //{data, allB, zjets, ttbar, lm4, zjets_un, wjets, znunu, vvjets, singletops, singletopt, singletopu, qcd100to250, qcd250to500, qcd500to1000, qcd1000toInf, BnoZ,BwithZ,BandS}
2315    
2316    
2317     for(int i=0;i<19;i++)
2318     {
2319     if(i==zjets_un)continue;
2320     if(i==BnoZ)continue;
2321     if(i==BwithZ)continue;
2322     if(i==BandS)continue;
2323    
2324     float countsError=0.;
2325     float counts = histIntegralAndError(ZPt_ossf[i],0,countsError);
2326    
2327    
2328     if(i==data)cout<<"data : ";
2329     if(i==allB)cout<<"allB : ";
2330     if(i==zjets)cout<<"zjets : ";
2331     if(i==ttbar)cout<<"ttbar : ";
2332     if(i==lm4)cout<<"lm4 : ";
2333     if(i==wjets)cout<<"wjets : ";
2334     if(i==znunu)cout<<"znunu : ";
2335     if(i==vvjets)cout<<"vvjets : ";
2336     if(i==singletops)cout<<"singletops : ";
2337     if(i==singletopt)cout<<"singletopt : ";
2338     if(i==singletopu)cout<<"singletopu : ";
2339     if(i==qcd100to250)cout<<"qcd100to250 : ";
2340     if(i==qcd250to500)cout<<"qcd250to500 : ";
2341     if(i==qcd500to1000)cout<<"qcd500to1000 : ";
2342     if(i==qcd1000toInf)cout<<"qcd1000toInf : ";
2343    
2344     cout<<setprecision(4)<<counts<<" +- "<<countsError<<endl;
2345     }
2346    
2347    
2348    
2349    
2350    
2351    
2352     }
2353    
2354    
2355    
2356     void setYMax(TH1F*hist,float x)
2357     {
2358     hist->SetMaximum(x*hist->GetMaximum());
2359     }
2360    
2361     void setYMin(TH1F*hist,float x)
2362     {
2363     hist->SetMinimum(x*hist->GetMinimum());
2364     }
2365    
2366     THStack *mc_stack(TH1F** histos)
2367     {
2368     THStack *hs_tmp = new THStack(histos[zjets]->GetName(),histos[zjets]->GetName());
2369     histos[znunu]->SetFillColor(kYellow);
2370     histos[qcd100to250]->SetFillColor(kPink);
2371     histos[singletops]->SetFillColor(kBlue);
2372     histos[singletopt]->SetFillColor(kBlue);
2373     histos[singletopu]->SetFillColor(kBlue);
2374     histos[wjets]->SetFillColor(kGray);
2375     histos[vvjets]->SetFillColor(kGreen);
2376     histos[ttbar]->SetFillColor(kMagenta+2);
2377     histos[zjets]->SetFillColor(kYellow);
2378    
2379     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2380     {
2381    
2382     int sampleIndex = allBList[counter];
2383     // histos[ttbar]->SetFillColor(kMagenta+2);
2384     // histos[zjets]->SetFillColor(kYellow);
2385     // if(sampleIndex!= ttbar && sampleIndex!=zjets)histos[sampleIndex]->SetFillColor(kBlack);
2386    
2387     // cout<<" sampleIndex = "<< sampleIndex<<endl;
2388     hs_tmp->Add(histos[sampleIndex]);
2389     }
2390    
2391     return hs_tmp;
2392     }
2393    
2394     THStack *mc_stack(TH1F** histos, int rebin) //blublu
2395     {
2396     stringstream stackname;
2397     stackname << histos[zjets]->GetName() << "__REBIN";
2398     THStack *hs_tmp = new THStack(stackname.str().c_str(),stackname.str().c_str());
2399     histos[znunu]->SetFillColor(kYellow);
2400     histos[qcd100to250]->SetFillColor(kPink);
2401     histos[singletops]->SetFillColor(kBlue);
2402     histos[singletopt]->SetFillColor(kBlue);
2403     histos[singletopu]->SetFillColor(kBlue);
2404     histos[wjets]->SetFillColor(kGray);
2405     histos[vvjets]->SetFillColor(kGreen);
2406     histos[ttbar]->SetFillColor(kMagenta+2);
2407     histos[zjets]->SetFillColor(kYellow);
2408    
2409     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2410     {
2411    
2412     int sampleIndex = allBList[counter];
2413     // histos[ttbar]->SetFillColor(kMagenta+2);
2414     // histos[zjets]->SetFillColor(kYellow);
2415     // if(sampleIndex!= ttbar && sampleIndex!=zjets)histos[sampleIndex]->SetFillColor(kBlack);
2416    
2417     // cout<<" sampleIndex = "<< sampleIndex<<endl;
2418     TH1F *reb = (TH1F*)histos[sampleIndex];
2419     reb->Rebin(rebin);
2420     hs_tmp->Add(reb);
2421     }
2422    
2423     return hs_tmp;
2424     }
2425    
2426     THStack *mc_normalizedtodata_stack(TH1F** histos, int rebin) //blublu
2427     {
2428     stringstream stackname;
2429     stackname << histos[zjets]->GetName() << "__REBIN";
2430     THStack *hs_tmp = new THStack(stackname.str().c_str(),stackname.str().c_str());
2431     histos[znunu]->SetFillColor(kYellow);
2432     histos[qcd100to250]->SetFillColor(kPink);
2433     histos[singletops]->SetFillColor(kBlue);
2434     histos[singletopt]->SetFillColor(kBlue);
2435     histos[singletopu]->SetFillColor(kBlue);
2436     histos[wjets]->SetFillColor(kGray);
2437     histos[vvjets]->SetFillColor(kGreen);
2438     histos[ttbar]->SetFillColor(kMagenta+2);
2439     histos[zjets]->SetFillColor(kYellow);
2440    
2441     float totalintegral=0;
2442     float dataintegral=histos[data]->Integral();
2443     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2444     {
2445    
2446     int sampleIndex = allBList[counter];
2447     float currintegral=histos[sampleIndex]->Integral();
2448     totalintegral+=currintegral;
2449     }
2450    
2451     cout << "------------------- " << histos[zjets]->GetName() << " Kostas Q2: What is the normalization? data:" << dataintegral << " vs mc int: " << totalintegral << " gives a ->Scale of " << dataintegral/totalintegral << endl;
2452     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2453     {
2454    
2455     int sampleIndex = allBList[counter];
2456     // histos[ttbar]->SetFillColor(kMagenta+2);
2457     // histos[zjets]->SetFillColor(kYellow);
2458     // if(sampleIndex!= ttbar && sampleIndex!=zjets)histos[sampleIndex]->SetFillColor(kBlack);
2459    
2460     // cout<<" sampleIndex = "<< sampleIndex<<endl;
2461     TH1F *reb = (TH1F*)histos[sampleIndex];
2462     reb->Rebin(rebin);
2463     reb->Scale(dataintegral/totalintegral);
2464     hs_tmp->Add(reb);
2465     }
2466    
2467     return hs_tmp;
2468     }
2469    
2470     THStack *mc_normalized_stack(TH1F** histos)
2471     {
2472     THStack *hs_tmp = new THStack(histos[zjets]->GetName(),histos[zjets]->GetName());
2473     histos[znunu]->SetFillColor(kYellow);
2474     histos[qcd100to250]->SetFillColor(kPink);
2475     histos[singletops]->SetFillColor(kBlue);
2476     histos[singletopt]->SetFillColor(kBlue);
2477     histos[singletopu]->SetFillColor(kBlue);
2478     histos[wjets]->SetFillColor(kGray);
2479     histos[vvjets]->SetFillColor(kGreen);
2480     histos[ttbar]->SetFillColor(kMagenta+2);
2481     histos[zjets]->SetFillColor(kYellow);
2482    
2483     float totalintegral=0;
2484     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2485     {
2486    
2487     int sampleIndex = allBList[counter];
2488     float currintegral=histos[sampleIndex]->Integral();
2489     totalintegral+=currintegral;
2490     }
2491     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2492     {
2493    
2494     int sampleIndex = allBList[counter];
2495     TH1F *currcop = (TH1F*)histos[sampleIndex]->Clone();
2496     currcop->Scale(1.0/totalintegral);
2497     //hs_tmp->Add(histos[sampleIndex]);
2498     hs_tmp->Add(currcop);
2499     }
2500    
2501     return hs_tmp;
2502     }
2503    
2504     THStack *mc_normalizedtodata_stack(TH1F** histos)
2505     {
2506    
2507     THStack *hs_tmp = new THStack(histos[zjets]->GetName(),histos[zjets]->GetName());
2508     histos[znunu]->SetFillColor(kYellow);
2509     histos[qcd100to250]->SetFillColor(kPink);
2510     histos[singletops]->SetFillColor(kBlue);
2511     histos[singletopt]->SetFillColor(kBlue);
2512     histos[singletopu]->SetFillColor(kBlue);
2513     histos[wjets]->SetFillColor(kGray);
2514     histos[vvjets]->SetFillColor(kGreen);
2515     histos[ttbar]->SetFillColor(kMagenta+2);
2516     histos[zjets]->SetFillColor(kYellow);
2517    
2518     float totalintegral=0;
2519     float dataintegral=histos[data]->Integral();
2520     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2521     {
2522    
2523     int sampleIndex = allBList[counter];
2524     float currintegral=histos[sampleIndex]->Integral();
2525     totalintegral+=currintegral;
2526     }
2527     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2528     {
2529    
2530     int sampleIndex = allBList[counter];
2531     TH1F *currcop = (TH1F*)histos[sampleIndex]->Clone();
2532     currcop->Scale(dataintegral/totalintegral);
2533     //hs_tmp->Add(histos[sampleIndex]);
2534     hs_tmp->Add(currcop);
2535     }
2536     cout << "------------------- " << histos[zjets]->GetName() << " Kostas Q2: What is the normalization? data:" << dataintegral << " vs mc int: " << totalintegral << " gives a ->Scale of " << dataintegral/totalintegral << endl;
2537    
2538     return hs_tmp;
2539     }
2540    
2541     TH1F *mc_th1f(TH1F** histos)
2542     {
2543     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2544     {
2545     int sampleIndex = allBList[counter];
2546     histos[allB]->Add(histos[sampleIndex]);
2547     }
2548    
2549     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2550     {
2551     int sampleIndex = allBList[counter];
2552     if(sampleIndex!=zjets)histos[BnoZ]->Add(histos[sampleIndex]);
2553     if(sampleIndex==zjets)histos[BwithZ]->Add(histos[sampleIndex]);
2554     }
2555    
2556     for(int counter=0 ; counter<int(allBList.size()) ; counter++)
2557     {
2558     int sampleIndex = allBList[counter];
2559     histos[BandS]->Add(histos[sampleIndex]);
2560     }
2561     histos[BandS]->Add(histos[lm4]);
2562    
2563     return histos[allB];
2564     }
2565    
2566    
2567     float histIntegral(TH1F *hist, float minMet)
2568     {
2569     // bin = 0; underflow bin
2570     // bin = 1; first bin with low-edge xlow INCLUDED
2571     // bin = nbins; last bin with upper-edge xup EXCLUDED
2572     // bin = nbins+1; overflow bin
2573    
2574     float val=0.;
2575     val = hist->Integral(hist->FindBin(double(minMet)),hist->GetNbinsX()+1);
2576     if(hist->GetBinContent(hist->GetNbinsX()+1)!=0)cout<< " hist " << hist->GetTitle() << " overflowed " << endl;
2577    
2578     return val;
2579     }
2580    
2581     float histIntegralAndError(TH1F *hist, float minMet, float &error)
2582     {
2583     // bin = 0; underflow bin
2584     // bin = 1; first bin with low-edge xlow INCLUDED
2585     // bin = nbins; last bin with upper-edge xup EXCLUDED
2586     // bin = nbins+1; overflow bin
2587    
2588     float val=0.;
2589     val = hist->Integral(hist->FindBin(double(minMet)),hist->GetNbinsX()+1);
2590     double valError=0.;
2591     hist->IntegralAndError(hist->FindBin(double(minMet)),hist->GetNbinsX()+1,valError);
2592     error = float(valError);
2593     return val;
2594     }
2595    
2596     float histIntegral(TH1F *hist, float minMet,float maxMet)
2597     {
2598     float val=0.;
2599     val = hist->Integral(hist->FindBin(double(minMet)),hist->FindBin(double(maxMet))+1);
2600     return val;
2601     }
2602    
2603     float histIntegralAndError(TH1F *hist, float minMet,float maxMet, float &error)
2604     {
2605     float val=0.;
2606     double valError=0.;
2607     val = hist->IntegralAndError(hist->FindBin(double(minMet)),hist->FindBin(double(maxMet))+1,valError);
2608     error = float(valError);
2609    
2610     // int bin1= hist->FindBin(double(minMet));
2611     // int bin2= hist->FindBin(double(maxMet))+1;
2612     // float sumw2=0;
2613     // for(int binx=bin1;binx<=bin2;++binx)
2614     // {
2615     // sumw2 += hist->GetBinError(binx)*hist->GetBinError(binx);
2616     // }
2617     // cout<<sqrt(sumw2)<<" error = "<<error<<endl;
2618    
2619     return val;
2620     }
2621    
2622     float computeRatioError(float a, float da, float b, float db)
2623     {
2624     float val=0.;
2625     float errorSquare = (a/b)*(a/b)*( (da/a)*(da/a) + (db/b)*(db/b));
2626     val = sqrt(errorSquare);
2627     return val;
2628    
2629     }
2630     float computeProductError(float a, float da, float b, float db)
2631     {
2632     float val=0.;
2633     float errorSquare = (a*b)*(a*b)*( (da/a)*(da/a) + (db/b)*(db/b));
2634     val = sqrt(errorSquare);
2635     return val;
2636     }
2637    
2638    
2639    
2640     template<class T>
2641     std::string any2string(T i)
2642     {
2643     std::ostringstream buffer;
2644     buffer << i;
2645     return buffer.str();
2646     }
2647    
2648    
2649     void newCanvas(int &canvasNumber)
2650     {
2651     can[canvasNumber] = new TCanvas(("c"+any2string(canvasNumber)).c_str(),("c"+any2string(canvasNumber)).c_str());
2652     can[canvasNumber]->cd();
2653     canvasNumber++;
2654     }
2655    
2656     void saveCanvas(int i)
2657     {
2658     string fileTitleC="fig/jzb_"+fileNamePreffix+"_"+any2string(i)+".C";
2659     string fileTitleE="fig/jzb_"+fileNamePreffix+"_"+any2string(i)+".eps";
2660     string fileTitleP="fig/jzb_"+fileNamePreffix+"_"+any2string(i)+".png";
2661     can[i]->SaveAs(fileTitleC.c_str());
2662     can[i]->SaveAs(fileTitleE.c_str());
2663     can[i]->SaveAs(fileTitleP.c_str());
2664     }
2665    
2666     void saveCanvas(int i, string filename)
2667     {
2668     string fileTitleC="fig/"+filename+".C";
2669     string fileTitleP="fig/"+filename+".png";
2670     string fileTitleE="fig/"+filename+".eps";
2671     can[i]->SaveAs(fileTitleC.c_str());
2672     can[i]->SaveAs(fileTitleP.c_str());
2673     can[i]->SaveAs(fileTitleE.c_str());
2674     }
2675    
2676    
2677    
2678     float get_peak(TH1F *hist, float &error, float &sigma, TF1* fitFunc, float lowlimit, float highlimit)
2679     {
2680     float mean = hist->GetBinCenter( hist->GetMaximumBin());
2681     float rms = hist->GetRMS();
2682    
2683     mean = hist->GetBinCenter( hist->GetMaximumBin());
2684    
2685     fitFunc->SetParameter(1,mean);
2686    
2687     hist->Fit(fitFunc,"QLL0","",mean-10,mean+10);
2688    
2689     mean=fitFunc->GetParameter(1);
2690     rms=fitFunc->GetParameter(2);
2691     error=fitFunc->GetParError(1);
2692    
2693     bool printOut = false; // print the peak estimate in the i-th iteration
2694    
2695     // --- perform iterations
2696     int numIterations=5;
2697    
2698     if(printOut) std::cout << " ( ";
2699     for(int i=1;i<numIterations+1;i++) //--- modify the number of iterations until peak is stable
2700     {
2701     hist->Fit(fitFunc,"QLLN","same",mean - lowlimit*rms, mean + highlimit*rms); // fit -2 +1 sigma from previous iteration
2702     mean=fitFunc->GetParameter(1);
2703     fitFunc->SetRange(mean - lowlimit*rms, mean + highlimit*rms);
2704     if(printOut) std::cout << mean << ",";
2705     }
2706     if(printOut) std::cout << " ) ";
2707     if(printOut) std::cout << endl;
2708     mean=fitFunc->GetParameter(1);
2709     sigma=fitFunc->GetParameter(2);
2710     error=1.0*fitFunc->GetParError(1);
2711    
2712     // cout << "[" << fitFunc->GetParameter(1) << " , " << fitFunc->GetParError(1) << "]" << endl;
2713     return mean;
2714     }
2715    
2716     TH1F *agreementHist(TH1F *ratioPlot)
2717     {
2718    
2719     string title = string(ratioPlot->GetTitle()) + "agreement";
2720     TH1F *agreement = new TH1F(title.c_str(),"pull distribution",100,0,3);
2721    
2722     int nBins = ratioPlot->GetNbinsX();
2723    
2724     for(int bin = 1; bin<=nBins+1 ; bin++)
2725     {
2726     float binContents = ratioPlot->GetBinContent(bin);
2727     agreement->Fill(binContents);
2728    
2729     }
2730    
2731     return agreement;
2732     }
2733    
2734     TH1F *agreementHist(TH1F *obs, TH1F *pred)
2735     {
2736    
2737     TH1F *agreement = (TH1F*)obs->Clone();
2738    
2739     int nBins = obs->GetNbinsX();
2740    
2741     for(int bin = 1; bin<=nBins+1 ; bin++)
2742     {
2743     float nobs = obs->GetBinContent(bin);
2744     float nobsError = obs->GetBinError(bin);
2745     float npred = pred->GetBinContent(bin);
2746     float npredError = pred->GetBinError(bin);
2747    
2748    
2749     float binContent = (nobs-npred)/sqrt(npredError*npredError + nobsError*nobsError);
2750     cout << nobs << " - " << npred << "binContent"<< binContent << endl;
2751     if(binContent<100)
2752     {
2753     agreement->SetBinContent(bin,binContent);
2754     }
2755     else
2756     {
2757     agreement->SetBinContent(bin,0);
2758     }
2759    
2760     }
2761    
2762     return agreement;
2763     }
2764    
2765     Double_t CrystalBall(Double_t *x,Double_t *par)
2766     {
2767     Double_t arg1=0,arg2=0,A=0,B=0;
2768     Double_t f1=0;
2769     Double_t f2=0;
2770     Double_t lim=0;
2771     Double_t fitval=0;
2772     Double_t N=0;
2773     Double_t n=par[4];
2774    
2775    
2776     if (par[2] != 0)
2777     arg1 = (x[0]-par[1])/par[2];
2778    
2779     arg2 = ( -pow( TMath::Abs(par[3]) , 2 ) ) / 2 ;
2780    
2781     if (par[3] != 0)
2782     A = pow( ( n / TMath::Abs( par[3] ) ) , n) * TMath::Exp(arg2);
2783    
2784     if (par[3] != 0)
2785     B = n / TMath::Abs(par[3]) - TMath::Abs(par[3]);
2786    
2787     f1 = TMath::Exp(-0.5*arg1*arg1);
2788     if (par[2] != 0)
2789     f2 = A * pow( ( B - (x[0] - par[1])/par[2] ) , -n );
2790    
2791     if (par[2] != 0)
2792     lim = ( par[1] - x[0] ) / par[2] ;
2793    
2794     N = par[0];
2795    
2796    
2797     if(lim < par[3])
2798     fitval = N * f1;
2799     if(lim >= par[3])
2800     fitval = N * f2;
2801    
2802     return fitval;
2803     }
2804    
2805    
2806    
2807     Double_t InvCrystalBall(Double_t *x,Double_t *par)
2808     {
2809     Double_t arg1=0,arg2=0,A=0,B=0;
2810     Double_t f1=0;
2811     Double_t f2=0;
2812     Double_t lim=0;
2813     Double_t fitval=0;
2814     Double_t N=0;
2815     Double_t n=par[4];
2816    
2817     Double_t invX = -x[0];
2818    
2819     if (par[2] != 0)
2820     arg1 = (invX-par[1])/par[2];
2821    
2822     arg2 = ( -pow( TMath::Abs(par[3]) , 2 ) ) / 2 ;
2823    
2824     if (par[3] != 0)
2825     A = pow( ( n / TMath::Abs( par[3] ) ) , n) * TMath::Exp(arg2);
2826    
2827     if (par[3] != 0)
2828     B = n / TMath::Abs(par[3]) - TMath::Abs(par[3]);
2829    
2830     f1 = TMath::Exp(-0.5*arg1*arg1);
2831     if (par[2] != 0)
2832     f2 = A * pow( ( B - (invX - par[1])/par[2] ) , -n );
2833    
2834     if (par[2] != 0)
2835     lim = ( par[1] - invX ) / par[2] ;
2836    
2837     N = par[0];
2838    
2839    
2840    
2841     if(lim < par[3])
2842     fitval = N * f1;
2843     if(lim >= par[3])
2844     fitval = N * f2;
2845    
2846     return fitval;
2847     }
2848    
2849     Double_t InvCrystalBallP(Double_t *x,Double_t *par)
2850     {
2851     Double_t arg1=0,arg2=0,A=0,B=0;
2852     Double_t f1=0;
2853     Double_t f2=0;
2854     Double_t lim=0;
2855     Double_t fitval=0;
2856     Double_t N=0;
2857     Double_t n=par[4];
2858    
2859     Double_t invX = -x[0];
2860    
2861     if (par[2] != 0)
2862     arg1 = (invX-par[1])/par[2];
2863    
2864     arg2 = ( -pow( TMath::Abs(par[3]) , 2 ) ) / 2 ;
2865    
2866     if (par[3] != 0)
2867     A = pow( ( n / TMath::Abs( par[3] ) ) , n) * TMath::Exp(arg2);
2868    
2869     if (par[3] != 0)
2870     B = n / TMath::Abs(par[3]) - TMath::Abs(par[3]);
2871    
2872     f1 = TMath::Exp(-0.5*arg1*arg1);
2873     if (par[2] != 0)
2874     f2 = A * pow( ( B - (invX - par[1])/par[2] ) , -n );
2875    
2876     if (par[2] != 0)
2877     lim = ( par[1] - invX ) / par[2] ;
2878    
2879     N = par[0];
2880    
2881    
2882    
2883     if(lim < par[3])
2884     fitval = N * f1;
2885     if(lim >= par[3])
2886     fitval = N * f2;
2887    
2888     fitval+= statErrorP(fitval);
2889     return fitval;
2890     }
2891    
2892     Double_t InvCrystalBallN(Double_t *x,Double_t *par)
2893     {
2894     Double_t arg1=0,arg2=0,A=0,B=0;
2895     Double_t f1=0;
2896     Double_t f2=0;
2897     Double_t lim=0;
2898     Double_t fitval=0;
2899     Double_t N=0;
2900     Double_t n=par[4];
2901    
2902     Double_t invX = -x[0];
2903    
2904     if (par[2] != 0)
2905     arg1 = (invX-par[1])/par[2];
2906    
2907     arg2 = ( -pow( TMath::Abs(par[3]) , 2 ) ) / 2 ;
2908    
2909     if (par[3] != 0)
2910     A = pow( ( n / TMath::Abs( par[3] ) ) , n) * TMath::Exp(arg2);
2911    
2912     if (par[3] != 0)
2913     B = n / TMath::Abs(par[3]) - TMath::Abs(par[3]);
2914    
2915     f1 = TMath::Exp(-0.5*arg1*arg1);
2916     if (par[2] != 0)
2917     f2 = A * pow( ( B - (invX - par[1])/par[2] ) , -n );
2918    
2919     if (par[2] != 0)
2920     lim = ( par[1] - invX ) / par[2] ;
2921    
2922     N = par[0];
2923    
2924    
2925    
2926     if(lim < par[3])
2927     fitval = N * f1;
2928     if(lim >= par[3])
2929     fitval = N * f2;
2930    
2931     fitval-= statErrorN(fitval);
2932     return fitval;
2933     }
2934    
2935     /*
2936     TH1F *histRatio(TH1F *h1,TH1F *h2, int id)
2937     {
2938     float absJZBbins[] = {0,5,10,20,35,55,80,110,200,480};
2939     int absJZBbinsNumber = sizeof(absJZBbins)/4-1;
2940    
2941     TH1F *tmpHist = new TH1F("JZBratio_allB","JZBratio_allB",absJZBbinsNumber,absJZBbins);
2942    
2943     for(int i=0;i<sizeof(absJZBbins)/4-1;i++)
2944     {
2945     float xmin = absJZBbins[i];
2946     float xmax = absJZBbins[i+1];
2947     float nominatorError = 0.;
2948     float nominator = histIntegralAndError(h1,xmin,xmax,nominatorError);
2949     float denominatorError=0.;
2950     float denominator = histIntegralAndError(h2,xmin,xmax,denominatorError);
2951     float error = computeRatioError(nominator,nominatorError,denominator,denominatorError);
2952     if(id==data)error = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
2953     if(denominator!=0)
2954     {
2955     tmpHist->SetBinContent(i+1,nominator/denominator);
2956     tmpHist->SetBinError(i+1,error);
2957     }
2958     }
2959    
2960    
2961     return tmpHist;
2962     }
2963    
2964     */
2965    
2966     TGraphAsymmErrors *histRatio(TH1F *h1,TH1F *h2, int id)
2967     {
2968     float absJZBbins[] = {0,5,10,20,35,50,75,100,350};
2969     // float absJZBbins[] = {0,5,10,20,35,55,80,110,200,480};
2970     // float absJZBbins[] = {0,5,10,20,30,40,50,60,70,80,90,110,200,480};
2971    
2972     int absJZBbinsNumber = sizeof(absJZBbins)/4-1;
2973     if(data)absJZBbinsNumber=absJZBbinsNumber-2; // do not bother for the last 2 bins, division by 0
2974    
2975     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
2976     // int nremovedbins = 0; // not in use anymore
2977    
2978     for(unsigned int i=0;i<sizeof(absJZBbins)/4-1;i++)
2979     {
2980     float xmin = absJZBbins[i];
2981     float xmax = absJZBbins[i+1];
2982     float xCenter = 0.5*(xmax+xmin);
2983     float xWidth = 0.5*(xmax-xmin);
2984     float nominatorError = 0.;
2985     float nominator = histIntegralAndError(h1,xmin,xmax-0.01,nominatorError);
2986     float denominatorError=0.;
2987     float denominator = histIntegralAndError(h2,xmin,xmax-0.01,denominatorError);
2988     float errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
2989     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
2990     if(id==data)
2991     {
2992     errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
2993     //errorN = computeRatioError(nominator,statErrorN(nominator),denominator,statErrorN(denominator));
2994     errorN = errorP; // symmetrize using statErrorP
2995     }
2996     if(denominator!=0)
2997     {
2998     cout << "[xmin,xmax] = ["<<xmin<<","<<xmax<<"] xCenter = "<<xCenter<<" xWidth = "<<xWidth<< " i = "<<i<<" value = "<<nominator/denominator<<endl;
2999     graph->SetPoint(i, xCenter, nominator/denominator);
3000     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
3001     }
3002     else
3003     {
3004     cout << "ARTIFICIALLY SET POINT: [xmin,xmax] = ["<<xmin<<","<<xmax<<"] xCenter = "<<xCenter<<" xWidth = "<<xWidth<< " i = "<<i<<" value = "<<nominator/denominator<<endl;
3005     graph->SetPoint(i, xCenter, nominator/denominator);
3006     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
3007     /*
3008     graph->RemovePoint(i-nremovedbins);
3009     nremovedbins++;
3010     cout<<"removing point"<<"[xmin,xmax] = ["<<xmin<<","<<xmax<<"] xCenter = "<<xCenter<<" xWidth = "<<xWidth<< " i = "<<i<<" denominator = "<<denominator<<" from id = "<< id <<endl;
3011     */
3012     }
3013     }
3014    
3015     return graph;
3016     }
3017    
3018    
3019     TGraphAsymmErrors *histRatiofin(TH1F *h1,TH1F *h2, int id)
3020     {
3021    
3022     int absJZBbinsNumber = sizeof(fin_binning)/4-1;
3023     // if(data)absJZBbinsNumber=absJZBbinsNumber-2; // do not bother for the last 2 bins, division by 0
3024    
3025     TGraphAsymmErrors* graph = new TGraphAsymmErrors(absJZBbinsNumber);
3026     // int nremovedbins = 0; // not in use anymore
3027    
3028     for(unsigned int i=0;i<sizeof(fin_binning)/4-1;i++)
3029     {
3030     /* float xmin = absJZBbins[i];
3031     float xmax = absJZBbins[i+1];
3032     float xCenter = 0.5*(xmax+xmin);
3033     float xWidth = 0.5*(xmax-xmin);*/
3034     // float xWidth=(fin_binning[i+1]-fin_binning[i])/2;
3035     // float xCenter=fin_binning[i]+0.5*xWidth;
3036     float xCenter=h1->GetBinCenter(i+1);
3037     float xWidth=(h1->GetBinWidth(i+1))*0.5;
3038     float nominatorError = h1->GetBinError(i+1);
3039     float nominator=h1->GetBinContent(i+1);
3040     // float nominator = histIntegralAndError(h1,xmin,xmax-0.01,nominatorError);
3041     float denominatorError=h2->GetBinError(i+1);
3042     float denominator=h2->GetBinContent(i+1);
3043     // float denominator = histIntegralAndError(h2,xmin,xmax-0.01,denominatorError);
3044     float errorN = computeRatioError(nominator,nominatorError,denominator,denominatorError);
3045     float errorP = computeRatioError(nominator,nominatorError,denominator,denominatorError);
3046     if(id==data)
3047     {
3048     errorP = computeRatioError(nominator,statErrorP(nominator),denominator,statErrorP(denominator));
3049     //errorN = computeRatioError(nominator,statErrorN(nominator),denominator,statErrorN(denominator));
3050     errorN = errorP; // symmetrize using statErrorP
3051     }
3052     if(denominator!=0)
3053     {
3054     // cout << "[xmin,xmax] = ["<<xmin<<","<<xmax<<"] xCenter = "<<xCenter<<" xWidth = "<<xWidth<< " i = "<<i<<" value = "<<nominator/denominator<<endl;
3055     graph->SetPoint(i, xCenter, nominator/denominator);
3056     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
3057     }
3058     else
3059     {
3060     // cout << "ARTIFICIALLY SET POINT: [xmin,xmax] = ["<<xmin<<","<<xmax<<"] xCenter = "<<xCenter<<" xWidth = "<<xWidth<< " i = "<<i<<" value = "<<nominator/denominator<<endl;
3061     graph->SetPoint(i, xCenter, nominator/denominator);
3062     graph->SetPointError(i,xWidth,xWidth,errorN,errorP);
3063     }
3064     }
3065    
3066     return graph;
3067     }
3068