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

# Content
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