source: huonglan/SPLOT/sPlots_NewNasty.C

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

re-upload sPlot version
File size: 10.9 KB
Line 
1double sP[12];
2//Set the cut values
3double Cut1 = 0.6;
4double Cut2 = 0.7;
5double Cut3 = 60;
6double Cut4 = 44;
7double Cut5 = 250;
8double Cut6 = 600;
9double sigmaMin = 0.167335;
10
11void sPlot(){
12  //######## COMPUTE the sWeights
13  TString filename = "";
14  TString type[2] = {"AllBkg","Signal"};
15  TFile *file[2];
16  TTree *t[2];
17  double WeiSig = 0.685*1.43*1000/74873;
18  TFile *fAllBkg = TFile::Open("res_AllBkg.root");
19  TTree *tAllBkg = (TTree*)fAllBkg->Get("T");
20  double WeiBkg = tAllBkg->GetWeight();
21  double Wei[2] = {WeiBkg,WeiSig};
22
23  //Define the tree for signal and background
24  for (int i = 0; i < 2; i++){
25    filename = "res_"; filename += type[i]; filename += ".root";
26    file[i] = TFile::Open(filename);
27    t[i] = (TTree*)file[i]->Get("T");
28  }
29 
30  double iCatS[12], iCatB[12];
31  double F1[12];
32  double F0[12];
33  int iTot[2];
34  for (int i = 0; i < 2; i++){
35    iTot[i] = 0;
36  }
37  for (int i = 0; i < 12; i++){
38    iCatS[i] = 0;
39    iCatB[i] = 0;
40  }
41
42  for (int i = 0; i < 2; i++){
43    int nevt = t[i]->GetEntries();
44    for (int iev = 0; iev < nevt; iev++) {
45      t[i]->GetEntry(iev);
46           
47      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
48      double HT    = t[i]->GetLeaf("HT")->GetValue();
49      double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
50      double y     = t[i]->GetLeaf("cosWqq")->GetValue();
51      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
52      double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
53     
54      double cos1 = fabs(x);
55      double cos2 = fabs(y);
56     
57      if (mt > Cut5 && HT > Cut6) {
58        iTot[i]++;
59
60        if (cos1<Cut1 && cos2<Cut2){
61          if ((dMLep<Cut3) && (dMHad<Cut4)){
62            if (i==0) iCatB[0]++;
63            if (i==1) iCatS[0]++;
64          }
65          if ((dMLep<Cut3) && (dMHad>Cut4)){
66            if (i==0) iCatB[1]++;
67            if (i==1) iCatS[1]++;
68          }
69          if ((dMLep>Cut3) && (dMHad<Cut4)){
70            if (i==0) iCatB[2]++;
71            if (i==1) iCatS[2]++;
72          }
73          if ((dMLep>Cut3) && (dMHad>Cut4)){
74            if (i==0) iCatB[3]++;
75            if (i==1) iCatS[3]++;
76          }
77        }
78
79        if (cos1>Cut1 && cos2>Cut2){
80          if ((dMLep<Cut3) && (dMHad<Cut4)){
81            if (i==0) iCatB[4]++;
82            if (i==1) iCatS[4]++;
83          }
84          if ((dMLep<Cut3) && (dMHad>Cut4)){
85            if (i==0) iCatB[5]++;
86            if (i==1) iCatS[5]++;
87          }
88          if ((dMLep>Cut3) && (dMHad<Cut4)){
89            if (i==0) iCatB[6]++;
90            if (i==1) iCatS[6]++;
91          }
92          if ((dMLep>Cut3) && (dMHad>Cut4)){
93            if (i==0) iCatB[7]++;
94            if (i==1) iCatS[7]++;
95          }
96        }
97       
98        if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
99          if ((dMLep<Cut3) && (dMHad<Cut4)){
100            if (i==0) iCatB[8]++;
101            if (i==1) iCatS[8]++;
102          }
103          if ((dMLep<Cut3) && (dMHad>Cut4)){
104            if (i==0) iCatB[9]++;
105            if (i==1) iCatS[9]++;
106          }
107          if ((dMLep>Cut3) && (dMHad<Cut4)){
108            if (i==0) iCatB[10]++;
109            if (i==1) iCatS[10]++;
110          }
111          if ((dMLep>Cut3) && (dMHad>Cut4)){
112            if (i==0) iCatB[11]++;
113            if (i==1) iCatS[11]++;
114          }
115        }
116      }
117    }//Read each file
118  }//End of reading 2 files
119 
120  double sumeffS = 0, sumeffB = 0;
121  for (int i = 0; i < 12; i++){
122    F1[i] = iCatS[i]/iTot[1];
123    F0[i] = iCatB[i]/iTot[0];
124    sumeffS += F1[i];
125    sumeffB += F0[i];
126  }
127
128  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
129
130  double N0 = iTot[0]*Wei[0];
131  double N1 = iTot[1]*Wei[1];
132
133  cout << "N1 = " << N1 << endl;
134  cout << "N0 = " << N0 << endl;
135
136  double sig2 = sigmaMin*sigmaMin*N1*N1;
137
138  for (int i = 0; i < 12; i++){
139    sP[i] = sig2*F1[i]/(N1*F1[i] + N0*F0[i]);
140  }
141}
142
143double sWeight(double cos1,double cos2,double dMLep,double dMHad,double Wei)
144{
145  double sWei = 0;
146
147  if (cos1<Cut1 && cos2<Cut2){
148    if ((dMLep<Cut3) && (dMHad<Cut4))
149      sWei = sP[0]*Wei;
150    if ((dMLep<Cut3) && (dMHad>Cut4))
151      sWei = sP[1]*Wei;
152    if ((dMLep>Cut3) && (dMHad<Cut4))
153      sWei = sP[2]*Wei;
154    if ((dMLep>Cut3) && (dMHad>Cut4))
155      sWei = sP[3]*Wei;
156  }
157 
158  if (cos1>Cut1 && cos2>Cut2){
159    if ((dMLep<Cut3) && (dMHad<Cut4))
160      sWei = sP[4]*Wei;
161    if ((dMLep<Cut3) && (dMHad>Cut4))
162      sWei = sP[5]*Wei;
163    if ((dMLep>Cut3) && (dMHad<Cut4))
164      sWei = sP[6]*Wei;
165    if ((dMLep>Cut3) && (dMHad>Cut4))
166      sWei = sP[7]*Wei;
167  }
168 
169  if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
170    if ((dMLep<Cut3) && (dMHad<Cut4))
171      sWei = sP[8]*Wei;
172    if ((dMLep<Cut3) && (dMHad>Cut4))
173      sWei = sP[9]*Wei;
174    if ((dMLep>Cut3) && (dMHad<Cut4))
175      sWei = sP[10]*Wei;
176    if ((dMLep>Cut3) && (dMHad>Cut4))
177      sWei = sP[11]*Wei;
178  }
179  return sWei;
180}
181
182//########################################################################
183//############################ MAIN PROGRAM ##############################
184//########################################################################
185
186void sPlots_NewNasty(){
187  gROOT->SetStyle("Plain");
188  //Make the table of sP values first
189  sPlot();
190
191  cout << "Cut1 = " << Cut1 << endl;
192  cout << "Cut2 = " << Cut2 << endl;
193  cout << "Cut3 = " << Cut3 << endl;
194  cout << "Cut4 = " << Cut4 << endl;
195  cout << "Cut5 = " << Cut5 << endl;
196  cout << "Cut6 = " << Cut6 << endl;
197
198  //Initialize the files and trees
199  TString filename = "";
200  TString type[2] = {"AllBkg","Signal"};
201  TFile *file[2];
202  TTree *t[2];
203  double WeiSig = 0.685*1.43*1000/74873;
204  TFile *fAllBkg = TFile::Open("res_AllBkg.root");
205  TTree *tAllBkg = (TTree*)fAllBkg->Get("T");
206  double WeiBkg = tAllBkg->GetWeight();
207  double Wei[2] = {WeiBkg,WeiSig};
208
209  //Define the tree for signal and background
210  for (int i = 0; i < 2; i++){
211    filename = "res_"; filename += type[i]; filename += ".root";
212    file[i] = TFile::Open(filename);
213    t[i] = (TTree*)file[i]->Get("T");
214  }
215 
216  TH2F *twoDim_mt_HT     = new TH2F("twoDim_mt_HT",     "", 150, 0, 1500, 100, 0, 2000); 
217  TH2F *twoDim_mt_HT_Bkg = new TH2F("twoDim_mt_HT_Bkg", "", 150, 0, 1500, 100, 0, 2000); 
218  TH2F *twoDim_mt_HT_Sig = new TH2F("twoDim_mt_HT_Sig", "", 150, 0, 1500, 100, 0, 2000); 
219
220  TH1F *hSig        = new TH1F("hSig",       "", 150, 0, 1500);
221  TH1F *hBkg        = new TH1F("hBkg",       "", 150, 0, 1500);
222  TH1F *h[10];
223  h[0] = new TH1F("hAll1ifb",  "",  150,  0,  1500);
224  h[1] = new TH1F("hSig1ifb",  "",  150,  0,  1500);
225  h[2] = new TH1F("hBkg1ifb",  "",  150,  0,  1500);
226  h[3] = new TH1F("sPlots",    "",  150,  0,  1500);
227  h[4] = new TH1F ("hSigW", "", 150, 0,  1500);
228  h[5] = new TH1F ("hBkgW", "", 150, 0,  1500);
229  h[6] = new TH1F("hBkgPure",  "",  150,  0,  1500);
230  h[7] = new TH1F("hSigFinal", "",  150,  0,  1500);
231  h[8] = new TH1F ("hsPlotsBkg", "", 150, 0,  1500);
232  h[9] = new TH1F ("hsPlotsSig", "", 150, 0,  1500);
233
234  //Compute the number of HyperNasty Events first
235  int iTot[2];
236  for (int i = 0; i < 2; i++){
237    iTot[i] = 0;
238  }
239  for (int i = 0; i < 2; i++){
240    int nevt = t[i]->GetEntries();
241    for (int iev = 0; iev < nevt; iev++) {
242      t[i]->GetEntry(iev);
243     
244      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
245      double HT   = t[i]->GetLeaf("HT")->GetValue();   
246      if (mt > Cut5 && HT > Cut6)  iTot[i]++;
247    }
248  }
249
250  double N1_bar = iTot[1]*Wei[1];
251  double N0_bar = iTot[0]*Wei[0];
252  double sig2 = sigmaMin*sigmaMin*N1_bar*N1_bar;
253
254  double k = (sig2 - N1_bar)/(N0_bar - (sig2 - N1_bar));
255
256  cout << "N1_bar = " << N1_bar << endl;
257  cout << "N0_bar = " << N0_bar << endl;
258
259  //start reading 2 files
260  for (int i = 0; i < 2; i++){
261    int nevt = t[i]->GetEntries();
262    for (int iev = 0; iev < nevt; iev++) {
263      t[i]->GetEntry(iev);
264     
265      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
266      double HT    = t[i]->GetLeaf("HT")->GetValue();
267      double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
268      double y     = t[i]->GetLeaf("cosWqq")->GetValue();
269      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
270      double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
271     
272      double cos1 = fabs(x);
273      double cos2 = fabs(y);
274
275      //sPlots total
276      if (mt > Cut5 && HT > Cut6){
277        double sWei  = sWeight(cos1,cos2,dMLep,dMHad,Wei[i]);
278        double sWeiN = sWei*(k+1) - k*Wei[i];
279        twoDim_mt_HT->Fill(mt,HT,sWei); 
280        if (i==0)
281          twoDim_mt_HT_Bkg->Fill(mt,HT,sWei);   
282
283        //1D plot
284        //For 1 fb-1
285        h[0]->Fill(mt,Wei[i]);
286        if (i==1) h[1]->Fill(mt,Wei[i]);
287        if (i==0) {
288          h[2]->Fill(mt,Wei[i]);
289          h[8]->Fill(mt,sWei);
290        }
291        //sPlots
292        h[3]->Fill(mt,sWei);
293        h[7]->Fill(mt,sWeiN);
294      } 
295    }//Read each file
296  }//End of reading 2 files
297
298  //twoDim_mt_HT_Bkg->Scale(1./twoDim_mt_HT_Bkg->Integral(0,1500,0,2000));
299  h[2]->Scale(1./h[2]->Integral(0,1500));
300
301  cout << "twoDim_mt_HT_Bkg integral is " << twoDim_mt_HT_Bkg->Integral(0,1500,0,2000) << endl;
302  cout << "h2 integral is " << h[2]->Integral(0,1500) << endl;
303
304  for (int i = 1; i < 151; i++){
305    double WSig = 0., WBkg = 0.;
306    for (int j = 1; j < 101; j++){
307      double vtot = twoDim_mt_HT->GetBinContent(i,j);
308      //double vBkg = (sig2 - N1_bar)*twoDim_mt_HT_Bkg->GetBinContent(i,j);
309      double vBkg = twoDim_mt_HT_Bkg->GetBinContent(i,j);
310
311      double vSig = vtot - vBkg;
312      twoDim_mt_HT_Sig->SetBinContent(i,j,vSig);
313      WSig += vSig;
314      WBkg += vBkg;
315    }
316    hSig->SetBinContent(i,WSig);
317    hBkg->SetBinContent(i,WBkg);
318  }
319
320  for (int i = 1; i < 151; i++){
321    double sP = h[3]->GetBinContent(i);
322    double B  = (sig2 - N1_bar)*h[2]->GetBinContent(i);
323    double sB = h[8]->GetBinContent(i);
324    double W  = sP - B;
325    double sW  = sP - sB;
326    double Tot1fb = h[0]->GetBinContent(i);
327    double BkgPure = Tot1fb - sP;
328
329    h[4]->SetBinContent(i,W);
330    h[9]->SetBinContent(i,sW);
331    h[5]->SetBinContent(i,B);
332    h[6]->SetBinContent(i,BkgPure);
333  }
334
335//   for (int i = 1; i < 151; i++){
336//     double v1 = hSig->GetBinContent(i);
337//     double v2 = h[4]->GetBinContent(i);
338//     double v3 = hBkg->GetBinContent(i);
339//     double v4 = h[5]->GetBinContent(i);
340//     cout << "v1 = " << v1 << " v2 = " << v2 << endl;
341//     cout << "v3 = " << v3 << " v4 = " << v4 << endl;
342//   }
343
344  cout << "h[4] integral is " << h[4]->Integral(0,1500) << endl;
345  cout << "h[5] integral is " << h[5]->Integral(0,1500) << endl;
346
347  TCanvas *c1 = new TCanvas("C1", "", 600, 500);
348  twoDim_mt_HT->Draw("lego2,0");
349
350  TCanvas *c2 = new TCanvas("C2", "", 600, 500);
351  twoDim_mt_HT_Bkg->Draw("lego2,0");
352
353  TCanvas *c3 = new TCanvas("C3", "", 600, 500);
354  twoDim_mt_HT_Sig->Draw("lego2,0");
355 
356  cout << "twoDim_mt_HT_Sig integral is " << twoDim_mt_HT_Sig->Integral(0,1500,0,2000) << endl; 
357
358  TCanvas *c4 = new TCanvas("C4", "", 1000, 500);
359  c4->Divide(4);
360  c4->cd(1);
361  h[4]->Draw();
362  c4->cd(2);
363  h[5]->Draw();
364  c4->cd(3);
365  h[9]->Draw();
366  c4->cd(4);
367  h[8]->Draw();
368  cout << "hSig integral is " << hSig->Integral(0,1500) << endl;
369  cout << "hBkg integral is " << hBkg->Integral(0,1500) << endl;
370  cout << "sig2 = " << sig2 << endl;
371  cout << "N1_bar = " << N1_bar << endl;
372  cout << "sig2 - N1_bar = " << sig2 - N1_bar << endl;
373
374  TCanvas *c5 = new TCanvas("C5", "", 1000, 800);
375  c5->Divide(3,2);
376  c5->cd(1);
377  h[0]->Draw();
378  c5->cd(2);
379  h[1]->Draw();
380  c5->cd(3);
381  h[2]->Draw();
382  c5->cd(4);
383  h[3]->Draw();
384  c5->cd(5);
385  h[4]->Draw();
386  c5->cd(6);
387  h[5]->Draw();
388
389  TCanvas *c6 = new TCanvas("C6", "", 800, 600);
390  c6->Divide(2);
391  c6->cd(1);
392  h[6]->Draw();
393  c6->cd(2);
394  h[7]->Draw();
395}
Note: See TracBrowser for help on using the repository browser.