source: huonglan/SPLOT/NEWREALEASE/sPlots_NewRel.C @ 19

Last change on this file since 19 was 19, checked in by huonglan, 13 years ago

sPlots for MC in the same file now, MakeFakeData with weight for trees contains scale factor for flavour fractions
File size: 18.7 KB
Line 
1#include <iostream>
2#include <TFile.h>
3#include <TTree.h>
4#include <TH1F.h>
5#include <TH2F.h>
6#include <TString.h>
7#include <TLeaf.h>
8#include <TRandom.h>
9#include <TMinuit.h>
10#include <TF1.h>
11#include <TCanvas.h>
12#include <vector>
13
14//Set the cut values
15double Cut1 = 0.7;
16double Cut2 = 0.7;
17double Cut3 = 48;
18double Cut4 = 250;
19double Cut5 = 620;
20double sigmaMin;
21int icount = 0;
22double iTot[2];
23double iCatS[6], iCatB[6];
24double F1[6];
25double F0[6];
26
27double N1 = 0., N0 = 0., sig2 = 0.;
28
29//Define class Event
30class EventC {
31public:
32  EventC() { m_w = 1; m_fs = 1; m_fb = 0; m_cat = 1; }
33  EventC(double w0, double fs0, double fb0) { m_w = w0; m_fs = fs0; m_fb = fb0; }
34  EventC(double w0, int cat0) { m_w = w0; m_cat = cat0; m_fs = 1; m_fb = 0; }
35  double w() { return m_w; }
36  double fs() { return m_fs; }
37  double fb() { return m_fb; }
38  int cat() { return m_cat; }
39  void fs(double fs0) { m_fs = fs0; }
40  void fb(double fb0) { m_fb = fb0; }
41
42private:
43  double m_w, m_fs, m_fb;
44  int m_cat;
45};
46
47std::vector<EventC> m_EventList;
48
49void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
50  double L = 0;
51  for (int iev = 0; iev < m_EventList.size(); iev++) {
52    EventC evt = m_EventList[iev];
53    double smt = par[0]*evt.fb() + par[1]*evt.fs();
54    if (smt > 0) L += evt.w()*log(smt);
55  }
56  f = par[0] + par[1] - L; 
57};
58
59
60vector<double> InvMatrix(double a, double b, double c, double d){
61  double detM = a*d - b*c;
62  double k = 1./detM;
63  double V00 = k*d;
64  double V01 = k*(-b);
65  double V10 = k*(-c);
66  double V11 = k*a;
67 
68  vector<double> Vinv;
69  Vinv.push_back(V00);
70  Vinv.push_back(V01);
71  Vinv.push_back(V10);
72  Vinv.push_back(V11);
73
74  return Vinv;
75}
76
77void CalF0F1(){
78  //######## COMPUTE the sWeights
79  TString filename = "";
80  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
81  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
82  TFile *file[2];
83  TTree *t[2];
84 
85  for (int i = 0; i < 2; i++){
86    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
87    file[i] = TFile::Open(filename);
88    t[i]    = (TTree*)file[i]->Get("T");
89  }
90
91  for (int i = 0; i < 2; i++){
92    iTot[i] = 0;
93  }
94  for (int i = 0; i < 6; i++){
95    iCatS[i] = 0;
96    iCatB[i] = 0;
97  }
98 
99  for (int i = 0; i < 2; i++){
100    int nevt = t[i]->GetEntries();
101    for (int iev = 0; iev < nevt; iev++) {
102      t[i]->GetEntry(iev);
103           
104      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
105      double HT     = t[i]->GetLeaf("HT")->GetValue();
106      double x      = t[i]->GetLeaf("cos1")->GetValue();
107      double y      = t[i]->GetLeaf("cos2")->GetValue();
108      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
109      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
110
111      //cout << "WeiTot = " << WeiTot << endl;
112
113      int hasBJ = t[i]->GetLeaf("hasBadJet")->GetValue();
114      double hf = t[i]->GetLeaf("hfor")->GetValue();
115
116      if (hf!=4 && hasBJ==0){
117       
118        double cos1 = fabs(x);
119        double cos2 = fabs(y);
120     
121        if (mt > Cut4 && HT > Cut5) {
122          icount++;
123          iTot[i] += WeiTot;
124          int it = -1;
125          if (cos1<Cut1 && cos2<Cut2){
126            if (dMLep<Cut3) it = 0;
127            if (dMLep>Cut3) it = 1;
128          }
129          else if (cos1>Cut1 && cos2>Cut2){
130            if (dMLep<Cut3) it = 2;
131            if (dMLep>Cut3) it = 3;
132          }
133          else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
134            if (dMLep<Cut3) it = 4;
135            if (dMLep>Cut3) it = 5;
136          } else {
137            cout << "Unknown category !!" << endl;
138          }
139          if (i == 0) 
140            iCatB[it] += WeiTot;
141          else
142            iCatS[it] += WeiTot;
143        }
144      }
145    }//Read each file
146  }//End of reading 2 files
147
148  cout << "iTot0 = " << iTot[0] << endl;
149  cout << "iTot1 = " << iTot[1] << endl;
150
151  N1 = iTot[1];
152  N0 = iTot[0];
153
154  double sigmaInv = 0;
155  double sumeffS = 0, sumeffB = 0;
156  for (int i = 0; i < 6; i++){
157    F1[i] = iCatS[i]/iTot[1];
158    F0[i] = iCatB[i]/iTot[0];
159    sumeffS += F1[i];
160    sumeffB += F0[i];
161    cout << "F1[" << i << "] = " << F1[i] << endl;
162    cout << "F0[" << i << "] = " << F0[i] << endl;
163
164    double term = F1[i]*F1[i]/(N1*F1[i] + N0*F0[i]);
165    if (F1[i]==0 && F0[i]==0)  term = 0;
166    sigmaInv += term;
167  }
168
169  double sigma2 = 1/N1/N1/sigmaInv;
170  sigmaMin = sqrt(sigma2);
171 
172  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
173  cout << "sigmaMin = " << sigmaMin << endl;
174}
175
176vector<double> vFe(double cos1,double cos2,double dMLep)
177{
178  double f0e = 0;
179  double f1e = 0;
180 
181  //cout << "cos1 = " << cos1 << " cos2 = " << cos2 << " dMLep = " << dMLep << endl;
182 
183  int it = -1;
184  if (cos1<Cut1 && cos2<Cut2){
185    if (dMLep<Cut3) it = 0;
186    if (dMLep>Cut3) it = 1;
187  }
188  else if (cos1>Cut1 && cos2>Cut2){
189    if (dMLep<Cut3) it = 2;
190    if (dMLep>Cut3) it = 3;
191  }
192  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
193    if (dMLep<Cut3) it = 4;
194    if (dMLep>Cut3) it = 5;
195  }
196 
197  f0e = F0[it];
198  f1e = F1[it];
199 
200  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
201 
202  vector<double> vfe;
203  vfe.push_back(f0e);
204  vfe.push_back(f1e);
205  return vfe;
206}
207
208void MakeEventListf0f1(){
209  //######## COMPUTE the sWeights
210  TString filename = "";
211  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
212  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
213  TFile *file[2];
214  TTree *t[2];
215 
216  for (int i = 0; i < 2; i++){
217    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
218    file[i] = TFile::Open(filename);
219    t[i]    = (TTree*)file[i]->Get("T");
220  }
221 
222  for (int i = 0; i < 2; i++){
223    int nevt = t[i]->GetEntries();
224    for (int iev = 0; iev < nevt; iev++) {
225      t[i]->GetEntry(iev);
226     
227      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
228      double HT     = t[i]->GetLeaf("HT")->GetValue();
229      double x      = t[i]->GetLeaf("cos1")->GetValue();
230      double y      = t[i]->GetLeaf("cos2")->GetValue();
231      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
232      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
233     
234      //cout << "WeiTot = " << WeiTot << endl;
235     
236      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
237      double hf = t[i]->GetLeaf("hfor")->GetValue();
238
239      if (hf!=4 && hasBJ==0){
240       
241        double cos1 = fabs(x);
242        double cos2 = fabs(y);
243       
244        if (mt > Cut4 && HT > Cut5) {
245          vector<double> fsfb =  vFe(cos1,cos2,dMLep);
246          double f0 = fsfb.at(0);
247          double f1 = fsfb.at(1);
248          //cout << "f0 = " << f0 << " f1 = " << f1 << endl;
249          EventC eventC(WeiTot,f1,f0);
250          m_EventList.push_back(eventC);
251        }
252      }//Read each file
253    }//End of reading 2 files
254  }
255}
256
257Fit(){
258  // Test
259  MakeEventListf0f1();
260
261  int ifl = 0;
262  double step = 0.01;
263  double fitPar0 = iTot[0];
264  double fitPar1 = iTot[1];
265  double nmax0   = fitPar0 + 4*sqrt(iTot[0]);
266  double nmin0   = fitPar0 - 4*sqrt(iTot[0]);
267  double nmax1   = fitPar1 + 4*sqrt(iTot[1]);
268  double nmin1   = - 4*sqrt(iTot[1]);
269  double n0best, err0best;
270  double n1best, err1best;
271  TMinuit* myMinuit = new TMinuit(2);
272  //myMinuit->Command("SET PRI -1");
273  myMinuit->Command("SET NOW");
274  myMinuit->Command("SET ERR 0.5");
275  myMinuit->Command("SET STR 2");   // try to improve mimimum search
276  myMinuit->SetFCN(fcn);
277  myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl); 
278  myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl); 
279  //myMinuit->FixParameter(0);
280  myMinuit->Command("MIGRAD 500 1");
281  myMinuit->GetParameter(0,n0best,err0best); 
282  myMinuit->GetParameter(1,n1best,err1best); 
283  std::cout << "Expected N0 = " << iTot[0] << " fitted = " << n0best << " +- " << err0best << std::endl;
284  std::cout << "Expected N1 = " << iTot[1] << " fitted = " << n1best << " +- " << err1best << std::endl;
285
286  N0 = n0best;
287  N1 = n1best;
288
289  cout << "N0 = " << N0 << " N1 = " << N1 << endl;
290  sig2 = sigmaMin*sigmaMin*N1*N1;
291}
292
293//########################################################################
294//############################ MAIN PROGRAM ##############################
295//########################################################################
296
297void sPlots_NewRel(){
298  //gROOT->SetStyle("Plain");
299  //Make the table of sP values first
300
301  cout << "Cut1 = " << Cut1 << endl;
302  cout << "Cut2 = " << Cut2 << endl;
303  cout << "Cut3 = " << Cut3 << endl;
304  cout << "Cut4 = " << Cut4 << endl;
305  cout << "Cut5 = " << Cut5 << endl;
306
307  CalF0F1();
308  cout << "icount = " << icount << endl;
309
310  //sPlots N0 known
311  cout << "***********************************************" << endl;
312  cout << "************** sPlot N0 known *****************" << endl;
313  cout << "***********************************************" << endl;
314  sPlots_NewRel_N0known();
315
316  //Get right N0 and N1 for standard sPlot
317  cout << "***********************************************" << endl;
318  cout << "************* sPlot N0 unknown ****************" << endl;
319  cout << "***********************************************" << endl;
320
321  Fit();
322 
323  TString filename = "";
324  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
325  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
326  TFile *file[2];
327  TTree *t[2];
328 
329  for (int i = 0; i < 2; i++){
330    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
331    file[i] = TFile::Open(filename);
332    t[i]    = (TTree*)file[i]->Get("T");
333  }
334
335  double VI00 = 0, VI01 = 0, VI10 = 0, VI11 = 0;
336  //start reading 2 files
337  //int icount = 0;
338  for (int i = 0; i < 2; i++){
339    int nevt = t[i]->GetEntries();
340    for (int iev = 0; iev < nevt; iev++) {
341      t[i]->GetEntry(iev);
342           
343      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
344      double HT     = t[i]->GetLeaf("HT")->GetValue();
345      double x      = t[i]->GetLeaf("cos1")->GetValue();
346      double y      = t[i]->GetLeaf("cos2")->GetValue();
347      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
348      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
349
350      double cos1 = fabs(x);
351      double cos2 = fabs(y);
352     
353      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
354      double hf = t[i]->GetLeaf("hfor")->GetValue();
355
356      if (hf!=4 && hasBJ==0){
357        if (mt > Cut4 && HT > Cut5){
358          //Compute the matrix
359          vector<double> veff = vFe(cos1,cos2,dMLep);
360          double F0e = veff.at(0);
361          double F1e = veff.at(1);
362         
363          //cout << "F0e = " << F0e << " F1e = " << F1e << endl;
364         
365          double vi00;
366          double vi01;
367          double vi10;
368          double vi11;
369         
370          vi00 = WeiTot*F0e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
371          vi01 = WeiTot*F0e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
372          vi10 = WeiTot*F1e*F0e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
373          vi11 = WeiTot*F1e*F1e/(N1*F1e + N0*F0e)/(N1*F1e + N0*F0e);
374
375          VI00 += vi00;
376          VI01 += vi01;
377          VI10 += vi10;
378          VI11 += vi11;
379        }
380      }
381    }//Read each file
382  }//End of reading 2 files
383
384  cout << "Main macro : N0 = " << N0 << " N1 = " << N1 << endl;
385
386  cout << "VI00 = " << VI00 << endl;
387  cout << "VI01 = " << VI01 << endl;
388  cout << "VI10 = " << VI10 << endl;
389  cout << "VI11 = " << VI11 << endl;
390  vector<double> Vinv = InvMatrix(VI00,VI01,VI10,VI11);
391  double V00 = Vinv.at(0);
392  double V01 = Vinv.at(1);
393  double V10 = Vinv.at(2);
394  double V11 = Vinv.at(3);
395  cout << "V00 = " << V00 << endl;
396  cout << "V01 = " << V01 << endl;
397  cout << "V10 = " << V10 << endl;
398  cout << "V11 = " << V11 << endl;
399  double sigmaN0 = sqrt(V00);
400  double sigmaN1 = sqrt(V11);
401  double rho     = V10/sqrt(V00*V11);
402  cout << "sigmaN0 = " << sigmaN0 << " sigmaN1 = " << sigmaN1 << " rho = " << rho << endl; 
403
404  TH1F *hSigN  = new TH1F("hSigN",     "", 35, 250, 600);
405  TH1F *hBkgN  = new TH1F("hBkgN",     "", 35, 250, 600);
406  TH2F *h2DSig = new TH2F("h2DSig",    "", 35, 250, 600, 50, 500, 1500);
407  TH2F *h2DBkg = new TH2F("h2DBkg",    "", 35, 250, 600, 50, 500, 1500);
408
409  //Try standard sPlots
410  double SumsP1 = 0.;
411  double SumsP0 = 0.;
412  double *sPN = new double[icount];
413  double *mtN = new double[icount];
414  int icc = 0;
415  for (int i = 0; i < 2; i++){
416    int nevt = t[i]->GetEntries();
417    for (int iev = 0; iev < nevt; iev++) {
418      t[i]->GetEntry(iev);
419     
420      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
421      double HT    = t[i]->GetLeaf("HT")->GetValue();
422      double x     = t[i]->GetLeaf("cos1")->GetValue();
423      double y     = t[i]->GetLeaf("cos2")->GetValue();
424      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
425      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
426     
427      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
428      double hf = t[i]->GetLeaf("hfor")->GetValue();
429
430      if (hf!=4 && hasBJ==0){
431       
432        double cos1 = fabs(x);
433        double cos2 = fabs(y);
434
435        //sPlots total
436        if (mt > Cut4 && HT > Cut5){
437          vector<double> veff = vFe(cos1,cos2,dMLep);
438          double F0e = veff.at(0);
439          double F1e = veff.at(1);
440         
441          double sP1 = (V11*F1e + V10*F0e)/(N1*F1e + N0*F0e);
442          double sP0 = (V00*F0e + V10*F1e)/(N1*F1e + N0*F0e);
443         
444          hSigN->Fill(mt,sP1*WeiTot);
445          hBkgN->Fill(mt,sP0*WeiTot);
446          h2DSig->Fill(mt,HT,sP1*WeiTot);
447          h2DBkg->Fill(mt,HT,sP0*WeiTot);
448          mtN[icc] = mt;
449          if (WeiTot==0) sPN[icc] = 0;
450          else sPN[icc] = sP1*sP1*WeiTot;
451          icc++;
452         
453          SumsP1 += sP1*WeiTot;
454          SumsP0 += sP0*WeiTot;
455        }
456      }
457    }//Read each file
458  }//End of reading 2 files
459 
460  cout << "SumsP0 = " << SumsP0 << " SumsP1 = " << SumsP1 << endl;
461
462  cout << endl << "*******************" << endl;
463  cout << "STANDARD SPLOTS" << endl;
464  cout << "hSigN integral  " << hSigN->Integral(0,1500)  << endl;
465  cout << "hBkgN integral  " << hBkgN->Integral(0,1500)  << endl;
466 
467  for (int i = 1; i < 36; i++){
468    double vb1 = 250 + (i-1)*10;
469    double vb2 = 250 + i*10;
470    double ErrBin = 0.;
471    for (int j = 0; j < icount; j++){
472      double mt = mtN[j];
473      double sp = sPN[j];
474      if (mt >= vb1 && mt < vb2){
475        ErrBin += sp;
476      }
477    }
478    hSigN->SetBinError(i,sqrt(ErrBin));
479  }
480
481  TH2F *frame = new TH2F("frame", "", 35, 250, 600, 50, -10, 40);
482
483  TCanvas *c = new TCanvas("c", "", 900, 400);
484  c->Divide(2);
485  c->cd(1); 
486  //frame->Draw();
487  hSigN->SetLineColor(2);
488  hSigN->SetLineWidth(1);
489  hSigN->DrawCopy("E1");
490  hSigN->SetLineColor(1);
491  hSigN->SetLineWidth(2.5);
492  hSigN->DrawCopy("hist,same");
493  c->cd(2);
494  hBkgN->SetLineColor(4);
495  hBkgN->SetLineWidth(2.5);
496  hBkgN->Draw();
497
498  TCanvas *c2D = new TCanvas("c2D", "", 900, 400);
499  c2D->Divide(2);
500  c2D->cd(1);
501  h2DSig->Draw("colz");
502  c2D->cd(2);
503  h2DBkg->Draw("colz");
504}
505
506
507//##################################################################
508//########### FUNCTION FOR BACKGROUND YIELD KNOWN ##################
509//##################################################################
510
511void sPlots_NewRel_N0known(){
512  cout << "in function N0 known " << icount << endl;
513  N0 = iTot[0];
514  N1 = iTot[1];
515  sig2 = sigmaMin*sigmaMin*N1*N1;
516  cout << "sigmaMin = " << sigmaMin << " sig2 = " << sig2 << endl;
517  cout << "N0 = " << N0 << " N1 = " << N1 << endl;
518
519  double k = (sig2 - N1)/(N0 - (sig2 - N1));
520 
521  TH2F *twoDim_mt_HT_Sig     = new TH2F("mt_HT_Sig",     "", 75, 250, 1500, 50, 500, 2000); 
522  TH2F *twoDim_mt_HT_Sig_esP = new TH2F("mt_HT_Sig_esP", "", 75, 250, 1500, 50, 500, 2000); 
523 
524  TH1F *h[7];
525  h[0]  = new TH1F("All1fb",       "",  150,  0,  1500);
526  h[1]  = new TH1F("Sig1fb",       "",  150,  0,  1500);
527  h[2]  = new TH1F("Bkg1fb",       "",  150,  0,  1500);
528  h[3]  = new TH1F("sPlotAll",     "",  150,  0,  1500);
529  h[4]  = new TH1F("sPlotSig",     "",  150,  0,  1500);
530  h[5]  = new TH1F("sPlotBkg",     "",  150,  0,  1500);
531  h[6]  = new TH1F("esPlotSig",    "",  150,  0,  1500);
532 
533  double *mtV  = new double[icount];
534  double *err  = new double[icount];
535  double *errN = new double[icount];
536  int icc = 0;
537 
538  TString filename = "";
539  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};
540  TFile *file[2];
541  TTree *t[2];
542  for (int i = 0; i < 2; i++){
543    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "NEW.root";
544    file[i] = TFile::Open(filename);
545    t[i]    = (TTree*)file[i]->Get("T");
546  }
547
548  for (int i = 0; i < 2; i++){
549    int nevt = t[i]->GetEntries();
550    for (int iev = 0; iev < nevt; iev++) {
551      t[i]->GetEntry(iev);
552     
553      double mt     = t[i]->GetLeaf("mtophad")->GetValue();
554      double HT     = t[i]->GetLeaf("HT")->GetValue();
555      double x      = t[i]->GetLeaf("cos1")->GetValue();
556      double y      = t[i]->GetLeaf("cos2")->GetValue();
557      double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
558      double WeiTot = t[i]->GetLeaf("wei")->GetValue();
559
560      int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
561      double hf = t[i]->GetLeaf("hfor")->GetValue();
562     
563      if (hf!=4 && hasBJ==0){
564       
565        double cos1 = fabs(x);
566        double cos2 = fabs(y);
567       
568        //sPlots total
569        if (mt > Cut4 && HT > Cut5){
570          vector<double> veff = vFe(cos1,cos2,dMLep);
571          double F0e = veff.at(0);
572          double F1e = veff.at(1);
573         
574          double sWei  = sig2*F1e*WeiTot/(N1*F1e + N0*F0e);
575          double sWeiN = sWei*(k+1) - k*WeiTot;
576         
577          //For 1 fb-1
578          h[0]->Fill(mt,WeiTot);
579          h[3]->Fill(mt,sWei);//sPlot
580          h[6]->Fill(mt,sWeiN);//esPlot
581
582          if (i==0) {
583            h[2]->Fill(mt,WeiTot);
584            h[5]->Fill(mt,sWei);
585          }
586          if (i==1){
587            h[1]->Fill(mt,WeiTot);
588            twoDim_mt_HT_Sig->Fill(mt,HT,sWei);
589            twoDim_mt_HT_Sig_esP->Fill(mt,HT,sWeiN);   
590          }
591         
592          //uncertainty
593          mtV[icc] = mt;
594          if (WeiTot==0){
595            err[icc] = 0;
596            errN[icc] = 0;
597          } else {
598            err[icc] = sWei*sWei/WeiTot;
599            errN[icc] = sWeiN*sWeiN/WeiTot;
600          }
601          icc++;
602        }
603      }
604    }//Read each file
605  }//End of reading 2 files
606 
607  h[2]->Scale(1./h[2]->Integral(0,1500));
608
609  for (int i = 1; i < 151; i++){
610    double sP = h[3]->GetBinContent(i);
611    double B  = (sig2 - N1)*h[2]->GetBinContent(i);
612    double sB = h[5]->GetBinContent(i);
613    double W  = sP - B;
614    double sW  = sP - sB;
615    double Tot1fb = h[0]->GetBinContent(i);
616    double BkgPure = Tot1fb - sP;
617   
618    double vb1 = 10*(i-1);
619    double vb2 = 10*i;
620    double SumErrBin = 0;
621    double SumErrNBin = 0;
622    for (int ic = 0; ic < icount; ic++){
623      double mt  = mtV[ic];
624      double errV  = err[ic];
625      double errNV = errN[ic];
626      if (mt>=vb1 && mt<vb2){
627        SumErrBin += errV;
628        SumErrNBin += errNV;
629      }
630    }
631    h[4]->SetBinContent(i,sW);
632    h[4]->SetBinError(i,sqrt(SumErrBin));
633    h[6]->SetBinError(i,sqrt(SumErrNBin));
634  }
635 
636  cout << "h2 integral is " << h[2]->Integral(0,1500) << endl;
637  cout << "h4 integral is " << h[4]->Integral(0,1000) << endl;
638  cout << "h6 integral is " << h[6]->Integral(0,1500) << endl;
639
640  double x1 = 250;
641  double x2 = 600 - 0.01;
642  for (int i = 0; i < 7; i++) {
643    //h[i]->Rebin(2);
644    h[i]->SetLineWidth(3);
645    TAxis *ax = h[i]->GetXaxis();
646    int b1 = ax->FindBin(x1);
647    int b2 = ax->FindBin(x2);
648    ax->SetRange(b1,b2);
649  }
650
651  //plot for sPlot with M0 known
652  TCanvas *csP = new TCanvas("csP", "", 800, 400);
653  csP->Divide(2);
654  csP->cd(1);
655  h[4]->SetLineColor(2);
656  h[4]->SetLineWidth(1);
657  h[4]->DrawCopy("E1");
658  h[1]->SetLineColor(kGreen+1);
659  h[1]->Draw("same");
660  h[4]->SetLineColor(1);
661  h[4]->SetLineWidth(3);
662  h[4]->DrawCopy("hist,same");
663  csP->cd(2);
664  twoDim_mt_HT_Sig->Draw("colz,0");
665
666  //plot for sPlot with M0 UNknown
667  TCanvas *cesP = new TCanvas("cesP", "", 800,400);
668  cesP->Divide(2);
669  cesP->cd(1);
670  h[6]->SetLineColor(2);
671  h[6]->SetLineWidth(1);
672  h[6]->DrawCopy("E1");
673  h[6]->SetLineColor(1);
674  h[6]->SetLineWidth(3);
675  h[6]->Draw("hist,same");
676  cesP->cd(2);
677  twoDim_mt_HT_Sig_esP->Draw("colz,0"); 
678}
Note: See TracBrowser for help on using the repository browser.