source: huonglan/SPLOT/NEWREALEASE/EFFICIENCY/EfficiencyDependOnCut.C @ 22

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

sPlots_NewRel improved
File size: 15.2 KB
Line 
1#include <TFile.h>
2#include <TH2F.h>
3#include <TTree.h>
4#include <TLeaf.h>
5#include <TChain.h>
6#include <TString.h>
7#include <TGraphErrors.h>
8#include <TGraph.h>
9#include <TCanvas.h>
10#include <TLegend.h>
11#include <TROOT.h>
12#include <TMinuit.h>
13
14#include <iostream>
15
16using namespace std;
17
18//Set the cut values
19double Cut1 = 0.7;
20double Cut2 = 0.7;
21double Cut3 = 48;
22double Cut4 = 250;
23double Cut5 = 620;
24
25//Define class Event
26class EventC {
27public:
28  EventC() { m_w = 1; m_fs = 1; m_fb = 0; m_cat = 1; }
29  EventC(double w0, double fs0, double fb0) { m_w = w0; m_fs = fs0; m_fb = fb0; }
30  EventC(double w0, int cat0) { m_w = w0; m_cat = cat0; m_fs = 1; m_fb = 0; }
31  double w() { return m_w; }
32  double fs() { return m_fs; }
33  double fb() { return m_fb; }
34  int cat() { return m_cat; }
35  void fs(double fs0) { m_fs = fs0; }
36  void fb(double fb0) { m_fb = fb0; }
37
38private:
39  double m_w, m_fs, m_fb;
40  int m_cat;
41};
42
43//global variables
44std::vector<EventC> m_EventList;
45double m_F1[3][31][6];
46double m_F0[3][31][6];
47double m_N1exp[3][31];
48double m_N0exp[3][31];
49
50void fcnN0N1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
51  double L = 0;
52  for (unsigned int iev = 0; iev < m_EventList.size(); iev++) {
53    EventC evt = m_EventList[iev];
54    double smt = par[0]*evt.fb() + par[1]*evt.fs();
55    if (smt > 0) L += evt.w()*log(smt);
56  }
57  f = par[0] + par[1] - L; 
58};
59
60vector<double> vFe(double cos1,double cos2,double dMLep, int i, int j)
61{
62  double f0e = 0;
63  double f1e = 0;
64 
65  //cout << "cos1 = " << cos1 << " cos2 = " << cos2 << " dMLep = " << dMLep << endl;
66 
67  int it = -1;
68  if (cos1<Cut1 && cos2<Cut2){
69    if (dMLep<Cut3) it = 0;
70    if (dMLep>Cut3) it = 1;
71  }
72  else if (cos1>Cut1 && cos2>Cut2){
73    if (dMLep<Cut3) it = 2;
74    if (dMLep>Cut3) it = 3;
75  }
76  else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
77    if (dMLep<Cut3) it = 4;
78    if (dMLep>Cut3) it = 5;
79  }
80 
81  f0e = m_F0[i][j][it];
82  f1e = m_F1[i][j][it];
83 
84  //cout << "it = " << it << " f0e = " << f0e << " f1e = " << f1e << endl;
85 
86  vector<double> vfe;
87  vfe.push_back(f0e);
88  vfe.push_back(f1e);
89  return vfe;
90}
91
92//Make list of events in Data file
93void MakeEventListf0f1(int i, int j){
94  m_EventList.clear(); 
95
96  TFile *fD = TFile::Open("/data/atlas/huonglan/FROMLYON/DATA/resData.root");
97  TTree *tD = (TTree*)fD->Get("T"); 
98  int nD = tD->GetEntries();
99  for (int iev = 0; iev < nD; iev++) {
100    tD->GetEntry(iev);
101   
102    double mt    = tD->GetLeaf("mtophad")->GetValue();
103    double HT    = tD->GetLeaf("HT")->GetValue();
104    double x     = tD->GetLeaf("cos_tlepCMtpair")->GetValue();
105    double y     = tD->GetLeaf("cosWqq")->GetValue();
106    double dMLep = tD->GetLeaf("MinLep1")->GetValue();
107
108    int hasBJ = (int)tD->GetLeaf("hasBadJet")->GetValue();
109    int RunN  = (int)tD->GetLeaf("RunNumber_ana")->GetValue();
110   
111    if (RunN >= 180614 && hasBJ==1) continue;
112
113    double cos1 = fabs(x);
114    double cos2 = fabs(y);
115   
116    //sPlots total
117    if (mt > Cut4 && HT > Cut5){
118      vector<double> fsfb =  vFe(cos1,cos2,dMLep,i,j);
119      double f0 = fsfb.at(0);
120      double f1 = fsfb.at(1);
121      //cout << "f0 = " << f0 << " f1 = " << f1 << endl;
122      EventC eventC(1,f1,f0);
123      m_EventList.push_back(eventC);
124    } 
125  }
126}
127
128//Fit the number N0 and N1 using the events from data files
129std::pair<pair<double,double>,double> Fit(double N0_exp, double N1_exp, int i, int j){
130  // Test
131  MakeEventListf0f1(i,j);
132  cout << "Number of events in the LH " << m_EventList.size() << endl;
133
134  int ifl = 0;
135  double step = 0.01;
136 
137  double fitPar0 = N0_exp;
138  double nmax0   = fitPar0 + 4*sqrt(N0_exp);
139  double nmin0   = fitPar0 - 4*sqrt(N0_exp);
140 
141  double fitPar1 = N1_exp;
142  double nmax1   = fitPar1 + 4*sqrt(N1_exp);
143  double nmin1   =         - 4*sqrt(N1_exp);
144 
145  double n0best, err0best;
146  double n1best, err1best;
147 
148  TMinuit* myMinuit = new TMinuit(2);
149  //myMinuit->Command("SET PRI -1");
150  myMinuit->Command("SET NOW");
151  myMinuit->Command("SET ERR 0.5");
152  myMinuit->Command("SET STR 2");   // try to improve mimimum search
153  myMinuit->SetFCN(fcnN0N1);
154  myMinuit->mnparm(0,"N0",fitPar0,step,nmin0,nmax0,ifl); 
155  myMinuit->mnparm(1,"N1",fitPar1,step,nmin1,nmax1,ifl); 
156  //myMinuit->FixParameter(0);
157  myMinuit->Command("MIGRAD 500 1");
158  myMinuit->GetParameter(0,n0best,err0best); 
159  myMinuit->GetParameter(1,n1best,err1best); 
160 
161  std::cout << "Expected N0 = " << N0_exp << " fitted = " << n0best << " +- " << err0best << std::endl;
162  std::cout << "Expected N1 = " << N1_exp << " fitted = " << n1best << " +- " << err1best << std::endl;
163 
164  double N0 = n0best;
165  double N1 = n1best;
166  double sig2 = err1best*err1best;
167
168  pair<double,double> N0N1 = make_pair<double,double> (N0,N1);
169 
170  return make_pair<pair<double,double>,double> (N0N1,sig2);
171};
172
173
174//##############################################################
175//####################### MAIN PROGRAM #########################
176//##############################################################
177
178void EfficiencyDependOnCut(){
179  TTree *t[4];
180
181  TFile *fB = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_AllBkg.root");
182  t[0] = (TTree*)fB->Get("T");
183
184  TFile *fS = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_Signal.root");
185  t[1] = (TTree*)fS->Get("T");
186
187  TFile *fBS = TFile::Open("../MERGEFILES/TESTEFFICIENCY/res_FakeData.root");
188  t[2] = (TTree*)fBS->Get("T");
189
190  TFile *fD = TFile::Open("/exp/atlas/huonglan/DATA2011/resData.root");
191  t[3] = (TTree*)fD->Get("T");
192
193  double   e1[3][4][31];
194  double   e2[3][4][31];   
195  double   e3[3][4][31];   
196  double   e4[3][4][31];   
197  double   e5[3][4][31];   
198  double   e6[3][4][31];   
199  double iTot[3][4][31];
200
201  for (int k = 0; k < 3; k++){
202    for (int i = 0; i < 4; i++){
203      for (int j = 0; j < 31; j++){
204        iTot[k][i][j] =0 ;
205        e1[k][i][j] = 0;
206        e2[k][i][j] = 0;
207        e3[k][i][j] = 0;
208        e4[k][i][j] = 0;
209        e5[k][i][j] = 0;
210        e6[k][i][j] = 0;
211      }
212    }
213  }
214
215  double mjjj[31];    for (int i = 0; i < 31; i++) { mjjj[i]    = 200 + i*10;}
216  double mjjjErr[31]; for (int i = 0; i < 31; i++) { mjjjErr[i] = 0;}
217
218  for (int i = 0; i < 4; i++){
219    int nevt = t[i]->GetEntries();
220    for (int iev = 0; iev < nevt; iev++) {
221      t[i]->GetEntry(iev);
222
223      //cout << "iev " << iev << endl;
224      double wei = 1;
225      double cos1 = 1, cos2 = 1;
226      int hfor = -1;
227
228      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
229      double ht   = t[i]->GetLeaf("HT")->GetValue();
230      int hasBJ   = (int)t[i]->GetLeaf("hasBadJet")->GetValue();
231      if (i==0||i==1||i==2) {
232        wei  = t[i]->GetLeaf("wei")->GetValue();
233        cos1 = fabs(t[i]->GetLeaf("cos1")->GetValue());
234        cos2 = fabs(t[i]->GetLeaf("cos2")->GetValue());
235        hfor = (int)t[i]->GetLeaf("hfor")->GetValue();
236      } else if (i==3){
237        cos1 = fabs(t[i]->GetLeaf("cos_tlepCMtpair")->GetValue());
238        cos2 = fabs(t[i]->GetLeaf("cosWqq")->GetValue());
239      }
240      double dM   = t[i]->GetLeaf("MinLep1")->GetValue();
241
242      if (hasBJ!=0 || hfor==4) continue;
243
244      int it = -1;
245      if (ht>620 && mt>200)           it = 0;
246      if (ht>500 && ht<620 && mt>200) it = 1;
247      if (ht>620 && mt>200 && mt<250) it = 2;
248
249      if (it==-1) continue;
250      //cout << "mt " << mt << " ht " << ht << " it " << it << endl;
251
252      for (int j = 0; j < 31; j++){
253        double mtCut = mjjj[j];
254
255        //cout << "mtCut " << mtCut << endl;
256         
257        if (mt>mtCut){
258          iTot[it][i][j] += wei;
259         
260          if (cos1<Cut1 && cos2<Cut2){
261            if (dM<Cut3) e1[it][i][j] += wei;
262            if (dM>Cut3) e2[it][i][j] += wei;
263          }
264          else if (cos1>Cut1 && cos2>Cut2){
265            if (dM<Cut3) e3[it][i][j] += wei;
266            if (dM>Cut3) e4[it][i][j] += wei;
267          }
268          else if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
269            if (dM<Cut3) e5[it][i][j] += wei;
270            if (dM>Cut3) e6[it][i][j] += wei;
271          }
272        }
273      }
274    }
275  }
276
277  double    eB[3][31][6],    eS[3][31][6],    eBS[3][31][6],    eD[3][31][6]; 
278  double eBErr[3][31][6], eSErr[3][31][6], eBSErr[3][31][6], eDErr[3][31][6]; 
279
280  for (int i = 0; i < 3; i++){
281    for (int j = 0; j < 31; j++){
282      //cout << "e1[" << i << "][0][" << j << "] " << e1[i][0][j];
283      //cout << " iTot[" << i << "][0][" << j << "] " << iTot[i][0][j] << endl;
284     
285      eB[i][j][0] = e1[i][0][j]/iTot[i][0][j];
286      eB[i][j][1] = e2[i][0][j]/iTot[i][0][j];
287      eB[i][j][2] = e3[i][0][j]/iTot[i][0][j];
288      eB[i][j][3] = e4[i][0][j]/iTot[i][0][j];
289      eB[i][j][4] = e5[i][0][j]/iTot[i][0][j];
290      eB[i][j][5] = e6[i][0][j]/iTot[i][0][j];
291
292      if (iTot[i][0][j]==0) { for (int k = 0; k < 6; k++) { eB[i][j][k] = 0; } }
293
294      eS[i][j][0] = e1[i][1][j]/iTot[i][1][j];
295      eS[i][j][1] = e2[i][1][j]/iTot[i][1][j];
296      eS[i][j][2] = e3[i][1][j]/iTot[i][1][j];
297      eS[i][j][3] = e4[i][1][j]/iTot[i][1][j];
298      eS[i][j][4] = e5[i][1][j]/iTot[i][1][j];
299      eS[i][j][5] = e6[i][1][j]/iTot[i][1][j];
300
301      if (iTot[i][1][j]==0) { for (int k = 0; k < 6; k++) { eS[i][j][k] = 0; } }
302     
303      eBS[i][j][0] = e1[i][2][j]/iTot[i][2][j];
304      eBS[i][j][1] = e2[i][2][j]/iTot[i][2][j];
305      eBS[i][j][2] = e3[i][2][j]/iTot[i][2][j];
306      eBS[i][j][3] = e4[i][2][j]/iTot[i][2][j];
307      eBS[i][j][4] = e5[i][2][j]/iTot[i][2][j];
308      eBS[i][j][5] = e6[i][2][j]/iTot[i][2][j];
309
310      if (iTot[i][2][j]==0) { for (int k = 0; k < 6; k++) { eBS[i][j][k] = 0; } }
311     
312      eD[i][j][0] = e1[i][3][j]/iTot[i][3][j]; 
313      eD[i][j][1] = e2[i][3][j]/iTot[i][3][j]; 
314      eD[i][j][2] = e3[i][3][j]/iTot[i][3][j]; 
315      eD[i][j][3] = e4[i][3][j]/iTot[i][3][j]; 
316      eD[i][j][4] = e5[i][3][j]/iTot[i][3][j]; 
317      eD[i][j][5] = e6[i][3][j]/iTot[i][3][j]; 
318
319      if (iTot[i][3][j]==0) { for (int k = 0; k < 6; k++) { eD[i][j][k] = 0; } }
320     
321      //cout << "e1S[" << i << "][" << j << "] = " << e1S[i][j];
322      //cout << " iTot[" << i << "][1][" << j << "] = " << iTot[i][1][j] << endl;
323    }
324  }
325 
326  //3 is the number of regions
327  //31 is the number of the cut for mjjj
328  //6 is the number of categories
329  for (int i = 0; i < 3; i++){
330    for (int j = 0; j < 31; j++){
331      for (int k = 0; k < 6; k++){
332        eBErr[i][j][k]  = sqrt(eB[i][j][k]*(1-eB[i][j][k])/iTot[i][0][j]);
333        eSErr[i][j][k]  = sqrt(eS[i][j][k]*(1-eS[i][j][k])/iTot[i][1][j]);
334        eBSErr[i][j][k] = sqrt(eBS[i][j][k]*(1-eBS[i][j][k])/iTot[i][2][j]);
335        eDErr[i][j][k]  = sqrt(eD[i][j][k]*(1-eD[i][j][k])/iTot[i][3][j]);
336
337        m_F1[i][j][k] = eS[i][j][k];
338        m_F0[i][j][k] = eB[i][j][k];
339        m_N1exp[i][j] = iTot[i][1][j];
340        m_N0exp[i][j] = iTot[i][0][j];
341
342        cout << "region " << i << " cut mjjj " << mjjj[j] << " cate " << k+1 ;
343        cout << " iTotB = " << iTot[i][0][j] ;
344        cout << " iTotS = " << iTot[i][1][j] << endl;
345      }
346    }
347  }
348
349  //3 is the number of regions
350  //31 is the number of the cut for mjjj
351  double N0fit[3][31], N1fit[3][31], N1err[3][31];
352  for (int i = 0; i < 3; i++){
353    for (int j = 0; j < 31; j++){
354      pair<pair<double,double>,double> N0N1sig2 = Fit(m_N0exp[i][j],m_N1exp[i][j],i,j);
355      pair<double,double> N0N1 = N0N1sig2.first;
356      double N0   = N0N1.first; 
357      double N1   = N0N1.second; 
358      double sig2 = N0N1sig2.second;
359      N0fit[i][j] = N0;
360      N1fit[i][j] = N1;
361      N1err[i][j] = sig2;
362    }
363  }
364
365  TH2F *frame1 =  new TH2F("Frame1", "", 31, 200, 500, 1000, 0.001, 1);
366  frame1->SetStats(0);
367  TH2F *frame2 =  new TH2F("Frame2", "", 31, 200, 500, 1000, 0.001, 1);
368  frame2->SetStats(0);
369  TH2F *frame3 =  new TH2F("Frame3", "", 31, 200, 500, 1000, 0.001, 1);
370  frame3->SetStats(0);
371
372  TLegend *legend = new TLegend(0.3,0.6,0.7,0.75);
373  legend->SetBorderSize(0);
374  legend->SetFillColor(0);
375
376  TString name = "";
377  TString region[3] = {"region1","region2","region3"};
378  TString cat[6] = {"cat1","cat2","cat3","cat4","cat5","cat6"};
379
380//   TGraphErrors *heffB[3][6];
381//   TGraphErrors *heffS[3][6];
382//   TGraphErrors *heffBS[3][6];
383//   TGraphErrors *heffD[3][6];
384
385  TGraph *heffB[3][6];
386  TGraph *heffS[3][6];
387  TGraph *heffBS[3][6];
388  TGraph *heffD[3][6];
389
390  //*************Define histograms on N1************
391  TGraphErrors *hN1[3];
392  for (int i = 0; i < 3; i++){
393    double N1fittemp[31], N1errtemp[31];
394    for (int j = 0; j < 31; j++){ 
395      N1fittemp[j] = N1fit[i][j];
396      N1errtemp[j] = N1err[i][j]; 
397      cout << "mjjj = " << mjjj[j] << endl;
398      cout << " m_N1exp[" << i << "][" << j << "]= " << m_N1exp[i][j] ;
399      cout << " m_N0exp[" << i << "][" << j << "]= " << m_N0exp[i][j] ;
400      cout << " N1fit[" << i << "][" << j << "]= " << N1fit[i][j] ;
401      cout << " N0fit[" << i << "][" << j << "]= " << N0fit[i][j] ;
402      cout << " N1err[" << i << "][" << j << "]= " << N1err[i][j] << endl;
403    }
404    hN1[i] = new TGraphErrors(31,mjjj,N1fittemp,mjjjErr,N1errtemp);
405  }
406
407  //*************Define histograms on efficiencies****************
408  for (int i = 0; i < 3; i++){
409    for (int j = 0; j < 6; j++){
410      double eff[31]; 
411      double effErr[31];
412
413      name = "heffB"; name += region[i]; name += cat[j];
414      for (int ic = 0; ic < 31; ic++) { 
415        eff[ic]    = eB[i][ic][j];
416        effErr[ic] = eBErr[i][ic][j];
417      }
418      //heffB[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
419      heffB[i][j] = new TGraph(31, mjjj, eff);
420      heffB[i][j]->SetLineColor(8);
421      heffB[i][j]->SetLineWidth(2);
422     
423      name = "heffS"; name += region[i]; name += cat[j];
424      for (int ic = 0; ic < 31; ic++) { 
425        eff[ic]    = eS[i][ic][j];
426        effErr[ic] = eSErr[i][ic][j];
427      }
428      //heffS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
429      heffS[i][j] = new TGraph(31, mjjj, eff);
430      heffS[i][j]->SetLineColor(2);
431      heffS[i][j]->SetLineWidth(2);
432     
433      name = "heffBS"; name += region[i]; name += cat[j];
434      for (int ic = 0; ic < 31; ic++) { 
435        eff[ic]    = eBS[i][ic][j];
436        effErr[ic] = eBSErr[i][ic][j];
437      }
438      //heffBS[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
439      heffBS[i][j] = new TGraph(31, mjjj, eff);
440      heffBS[i][j]->SetLineColor(kMagenta);
441      heffBS[i][j]->SetLineWidth(2);
442     
443      name = "heffD"; name += region[i]; name += cat[j];
444      for (int ic = 0; ic < 31; ic++) { 
445        eff[ic]    = eD[i][ic][j];
446        effErr[ic] = eDErr[i][ic][j];
447      }
448      //heffD[i][j] = new TGraphErrors(31, mjjj, eff, mjjjErr, effErr);
449      heffD[i][j] = new TGraph(31, mjjj, eff);
450      heffD[i][j]->SetLineWidth(2);
451    }
452  }
453
454  legend->AddEntry(heffB[0][1],"BkgMC");
455  legend->AddEntry(heffS[0][1],"SigMC");
456  legend->AddEntry(heffBS[0][1],"Sig+BkgMC");
457  legend->AddEntry(heffD[0][1],"Data");
458
459  //****************** DRAWING *********************
460
461  TCanvas *cN1 = new TCanvas("cN1", "cN1", 1000, 400);
462  cN1->Divide(3);
463 
464  cN1->cd(1);
465  hN1[0]->Draw("ALP");
466  cN1->cd(2);
467  hN1[1]->Draw("ALP");
468  cN1->cd(3);
469  hN1[2]->Draw("ALP");
470
471  cN1->Print("N1fit.eps");
472
473  TCanvas *c[6];
474  for (int i = 0; i < 6; i++){
475    name = "Cat"; name += (i+1);
476    c[i] = new TCanvas(name, name, 1000, 400);
477    c[i]->Divide(3);
478
479    c[i]->cd(1);
480    name += ": mt>250 && HT>620";
481    frame1->SetTitle(name);
482    frame1->Draw();
483    heffS[0][i]->Draw("LP");
484    heffD[0][i]->Draw();
485    heffBS[0][i]->Draw();
486    heffB[0][i]->Draw();
487    legend->Draw("same");
488
489    c[i]->cd(2);
490    name = "mt>200 && HT>500 && HT<620";
491    frame2->SetTitle(name);
492    frame2->Draw();
493    heffS[1][i]->Draw("LP");
494    heffD[1][i]->Draw();
495    heffBS[1][i]->Draw();
496    heffB[1][i]->Draw();
497    legend->Draw("same");
498
499    c[i]->cd(3);
500    name = "mt>200 && mt<250 && HT>620";
501    frame3->SetTitle(name);
502    frame3->Draw();
503    heffS[2][i]->Draw("LP");
504    heffD[2][i]->Draw();
505    heffBS[2][i]->Draw();
506    heffB[2][i]->Draw();
507    legend->Draw("same");
508
509    name = "Cat"; name += (i+1); name += ".eps";
510    c[i]->Print(name);
511  }
512
513}
Note: See TracBrowser for help on using the repository browser.