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

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

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