source: huonglan/SPLOT/sPlots_2DeltaMin_NoHT.C

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

re-upload sPlot version
File size: 8.1 KB
Line 
1double sP[12];
2//Set the cut values
3double Cut1 = 0.6;
4double Cut2 = 0.7;
5double Cut3 = 54;
6double Cut4 = 44;
7double sigmaMin = 1.76824;
8
9void sPlot(){
10  //######## COMPUTE the sWeights
11  TString filename = "";
12  TString type[2] = {"AllBkg","Signal"};
13  TFile *file[2];
14  TTree *t[2];
15  double WeiSig = 0.685*1.43*1000/74873;
16  TFile *fAllBkg = TFile::Open("res_AllBkg.root");
17  TTree *tAllBkg = (TTree*)fAllBkg->Get("T");
18  double WeiBkg = tAllBkg->GetWeight();
19  double Wei[2] = {WeiBkg,WeiSig};
20
21  //Define the tree for signal and background
22  for (int i = 0; i < 2; i++){
23    filename = "res_"; filename += type[i]; filename += ".root";
24    file[i] = TFile::Open(filename);
25    t[i] = (TTree*)file[i]->Get("T");
26  }
27 
28  double iCatS[12], iCatB[12];
29  double F1[12];
30  double F0[12];
31  int iTot[2];
32  for (int i = 0; i < 2; i++){
33    iTot[i] = 0;
34  }
35  for (int i = 0; i < 12; i++){
36    iCatS[i] = 0;
37    iCatB[i] = 0;
38  }
39
40  for (int i = 0; i < 2; i++){
41    int nevt = t[i]->GetEntries();
42    for (int iev = 0; iev < nevt; iev++) {
43      t[i]->GetEntry(iev);
44           
45      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
46      double HT    = t[i]->GetLeaf("HT")->GetValue();
47      double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
48      double y     = t[i]->GetLeaf("cosWqq")->GetValue();
49      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
50      double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
51     
52      double cos1 = fabs(x);
53      double cos2 = fabs(y);
54     
55      if (mt > 300 && HT > 600) {
56        iTot[i]++;
57
58        if (cos1<Cut1 && cos2<Cut2){
59          if ((dMLep<Cut3) && (dMHad<Cut4)){
60            if (i==0) iCatB[0]++;
61            if (i==1) iCatS[0]++;
62          }
63          if ((dMLep<Cut3) && (dMHad>Cut4)){
64            if (i==0) iCatB[1]++;
65            if (i==1) iCatS[1]++;
66          }
67          if ((dMLep>Cut3) && (dMHad<Cut4)){
68            if (i==0) iCatB[2]++;
69            if (i==1) iCatS[2]++;
70          }
71          if ((dMLep>Cut3) && (dMHad>Cut4)){
72            if (i==0) iCatB[3]++;
73            if (i==1) iCatS[3]++;
74          }
75        }
76
77        if (cos1>Cut1 && cos2>Cut2){
78          if ((dMLep<Cut3) && (dMHad<Cut4)){
79            if (i==0) iCatB[4]++;
80            if (i==1) iCatS[4]++;
81          }
82          if ((dMLep<Cut3) && (dMHad>Cut4)){
83            if (i==0) iCatB[5]++;
84            if (i==1) iCatS[5]++;
85          }
86          if ((dMLep>Cut3) && (dMHad<Cut4)){
87            if (i==0) iCatB[6]++;
88            if (i==1) iCatS[6]++;
89          }
90          if ((dMLep>Cut3) && (dMHad>Cut4)){
91            if (i==0) iCatB[7]++;
92            if (i==1) iCatS[7]++;
93          }
94        }
95       
96        if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
97          if ((dMLep<Cut3) && (dMHad<Cut4)){
98            if (i==0) iCatB[8]++;
99            if (i==1) iCatS[8]++;
100          }
101          if ((dMLep<Cut3) && (dMHad>Cut4)){
102            if (i==0) iCatB[9]++;
103            if (i==1) iCatS[9]++;
104          }
105          if ((dMLep>Cut3) && (dMHad<Cut4)){
106            if (i==0) iCatB[10]++;
107            if (i==1) iCatS[10]++;
108          }
109          if ((dMLep>Cut3) && (dMHad>Cut4)){
110            if (i==0) iCatB[11]++;
111            if (i==1) iCatS[11]++;
112          }
113        }
114      }
115    }//Read each file
116  }//End of reading 2 files
117 
118  double sumeffS = 0, sumeffB = 0;
119  for (int i = 0; i < 12; i++){
120    F1[i] = iCatS[i]/iTot[1];
121    F0[i] = iCatB[i]/iTot[0];
122    sumeffS += F1[i];
123    sumeffB += F0[i];
124  }
125
126  cout << "sumeffS = " << sumeffS << " sumeffB = " << sumeffB << endl;
127
128  double N0 = iTot[0]*Wei[0];
129  double N1 = iTot[1]*Wei[1];
130
131  cout << "N1 = " << N1 << endl;
132  cout << "N0 = " << N0 << endl;
133
134  double sig2 = sigmaMin*sigmaMin*N1;
135
136  for (int i = 0; i < 12; i++){
137    sP[i] = sig2*F1[i]/(N1*F1[i] + N0*F0[i]);
138  }
139}
140
141double sWeight(double cos1,double cos2,double dMLep,double dMHad,double Wei)
142{
143  double sWei = 0;
144
145  if (cos1<Cut1 && cos2<Cut2){
146    if ((dMLep<Cut3) && (dMHad<Cut4))
147      sWei = sP[0]*Wei;
148    if ((dMLep<Cut3) && (dMHad>Cut4))
149      sWei = sP[1]*Wei;
150    if ((dMLep>Cut3) && (dMHad<Cut4))
151      sWei = sP[2]*Wei;
152    if ((dMLep>Cut3) && (dMHad>Cut4))
153      sWei = sP[3]*Wei;
154  }
155 
156  if (cos1>Cut1 && cos2>Cut2){
157    if ((dMLep<Cut3) && (dMHad<Cut4))
158      sWei = sP[4]*Wei;
159    if ((dMLep<Cut3) && (dMHad>Cut4))
160      sWei = sP[5]*Wei;
161    if ((dMLep>Cut3) && (dMHad<Cut4))
162      sWei = sP[6]*Wei;
163    if ((dMLep>Cut3) && (dMHad>Cut4))
164      sWei = sP[7]*Wei;
165  }
166 
167  if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
168    if ((dMLep<Cut3) && (dMHad<Cut4))
169      sWei = sP[8]*Wei;
170    if ((dMLep<Cut3) && (dMHad>Cut4))
171      sWei = sP[9]*Wei;
172    if ((dMLep>Cut3) && (dMHad<Cut4))
173      sWei = sP[10]*Wei;
174    if ((dMLep>Cut3) && (dMHad>Cut4))
175      sWei = sP[11]*Wei;
176  }
177  return sWei;
178}
179
180//########################################################################
181//############################ MAIN PROGRAM ##############################
182//########################################################################
183
184void sPlots_2DeltaMin_NoHT(){
185  gROOT->SetStyle("Plain");
186  //Make the table of sP values first
187  sPlot();
188
189  cout << "Cut1 = " << Cut1 << endl;
190  cout << "Cut2 = " << Cut2 << endl;
191  cout << "Cut3 = " << Cut3 << endl;
192  cout << "Cut4 = " << Cut4 << endl;
193
194  //Initialize the files and trees
195  TString filename = "";
196  TString type[2] = {"AllBkg","Signal"};
197  TFile *file[2];
198  TTree *t[2];
199  double WeiSig = 0.685*1.43*1000/74873;
200  TFile *fAllBkg = TFile::Open("res_AllBkg.root");
201  TTree *tAllBkg = (TTree*)fAllBkg->Get("T");
202  double WeiBkg = tAllBkg->GetWeight();
203  double Wei[2] = {WeiBkg,WeiSig};
204
205  //Define the tree for signal and background
206  for (int i = 0; i < 2; i++){
207    filename = "res_"; filename += type[i]; filename += ".root";
208    file[i] = TFile::Open(filename);
209    t[i] = (TTree*)file[i]->Get("T");
210  }
211 
212  TH2F *twoDim_mt_HT     = new TH2F("twoDim_mt_HT",     "", 150, 0, 1500, 100, 0, 2000); 
213  TH2F *twoDim_mt_HT_Bkg = new TH2F("twoDim_mt_HT_Bkg", "", 150, 0, 1500, 100, 0, 2000); 
214  TH2F *twoDim_mt_HT_Sig = new TH2F("twoDim_mt_HT_Sig", "", 150, 0, 1500, 100, 0, 2000); 
215
216  TH1F *hSig        = new TH1F("hSig",       "", 150, 0, 1500);
217
218  //Compute the number of HyperNasty Events first
219  int iTot[2];
220  for (int i = 0; i < 2; i++){
221    iTot[i] = 0;
222  }
223  for (int i = 0; i < 2; i++){
224    int nevt = t[i]->GetEntries();
225    for (int iev = 0; iev < nevt; iev++) {
226      t[i]->GetEntry(iev);
227     
228      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
229      double HT   = t[i]->GetLeaf("HT")->GetValue();   
230      if (mt > 300 && HT > 600)  iTot[i]++;
231    }
232  }
233
234  double N1_bar = iTot[1]*Wei[1];
235  double N0_bar = iTot[0]*Wei[0];
236
237  cout << "N1_bar = " << N1_bar << endl;
238  cout << "N0_bar = " << N0_bar << endl;
239
240  double sig2 = sigmaMin*sigmaMin*N1_bar;
241
242  //start reading 2 files
243  for (int i = 0; i < 2; i++){
244    int nevt = t[i]->GetEntries();
245    for (int iev = 0; iev < nevt; iev++) {
246      t[i]->GetEntry(iev);
247     
248      double mt    = t[i]->GetLeaf("mtophad")->GetValue();
249      double HT    = t[i]->GetLeaf("HT")->GetValue();
250      double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
251      double y     = t[i]->GetLeaf("cosWqq")->GetValue();
252      double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
253      double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
254     
255      double cos1 = fabs(x);
256      double cos2 = fabs(y);
257
258      //sPlots total
259      if (mt>300&&HT>600){
260        double sWei  = sWeight(cos1,cos2,dMLep,dMHad,Wei[i]);
261        twoDim_mt_HT->Fill(mt,HT,sWei); 
262        if (i==0)
263          twoDim_mt_HT_Bkg->Fill(mt,HT,sWei);   
264      } 
265    }//Read each file
266  }//End of reading 2 files
267
268  twoDim_mt_HT_Bkg->Scale(1./twoDim_mt_HT_Bkg->Integral(0,1500,0,2000));
269
270  cout << "twoDim_mt_HT_Bkg integral is " << twoDim_mt_HT_Bkg->Integral(0,1500,0,2000) << endl;
271
272  for (int i = 0; i < 150; i++){
273    double Wmt = 0.;
274    for (int j = 0; j < 100; j++){
275      double vtot = twoDim_mt_HT->GetBinContent(i,j);
276      double vBkg = twoDim_mt_HT_Bkg->GetBinContent(i,j);
277
278      double vSig = vtot - (sig2 - N1_bar)*vBkg;
279      twoDim_mt_HT_Sig->SetBinContent(i,j,vSig);
280      Wmt += vSig;
281    }
282    hSig->SetBinContent(i,Wmt);
283  }
284
285  TCanvas *c1 = new TCanvas("C1", "", 600, 500);
286  twoDim_mt_HT->Draw("lego2,0");
287
288  TCanvas *c2 = new TCanvas("C2", "", 600, 500);
289  twoDim_mt_HT_Bkg->Draw("lego2,0");
290
291  TCanvas *c3 = new TCanvas("C3", "", 600, 500);
292  twoDim_mt_HT_Sig->Draw("lego2,0");
293 
294  cout << "twoDim_mt_HT_Sig integral is " << twoDim_mt_HT_Sig->Integral(0,1500,0,2000) << endl; 
295
296  TCanvas *c4 = new TCanvas("C4", "", 600, 500);
297  hSig->Draw();
298  cout << "hSig integral is " << hSig->Integral(0,150) << endl;
299}
Note: See TracBrowser for help on using the repository browser.