source: huonglan/SPLOT/OptimizeCutWithoutMinLep.C

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

re-upload sPlot version
File size: 7.4 KB
Line 
1{
2  gROOT->SetStyle("Plain");
3  TString filename = "";
4  TString type[2] = {"Signal","ttbar"};
5  TFile *file[2];
6  TTree *t[2];
7  double WeiSig = 0.685*1.43*1000/74873;
8  double WeiTBkg = 0.5558*164.6*1000/199869;
9  double Wei[2] = {WeiSig,WeiTBkg};
10 
11  FILE *fRes3 = fopen("CatEfficiency_WithoutMinLep.txt","w+");
12 
13  for (int i = 0; i < 2; i++){
14    filename = "res_medium2ndEl_NewMinLep1_"; filename += type[i]; filename += ".root";
15    file[i] = TFile::Open(filename);
16    t[i] = (TTree*)file[i]->Get("T");
17    //t[i]->SetWeight(Wei[i],"global");
18  }
19
20  //#########   OPTIMIZE THE CUTS   ##########
21  double cos1V[11];
22  double cos2V[11];
23  double dMinV[1];
24  double C1Best = -1e3;
25  double C2Best = -1e3;
26  double C3Best = -1e3;
27  double Rmax   = -1e30;
28  double nSig   = -1e3;
29  double nBkg   = -1e3;
30
31  int iTot[2];
32  for (int i = 0; i < 2; i++){
33    iTot[i] = 0;
34  }
35
36  int i1Sig[121];
37  int i2Sig[121];
38  int i3Sig[121];
39  int i4Sig[121];
40  int i5Sig[121];
41  int i6Sig[121];
42
43  int i1Bkg[121];
44  int i2Bkg[121];
45  int i3Bkg[121];
46  int i4Bkg[121];
47  int i5Bkg[121];
48  int i6Bkg[121];
49
50  for (int i = 0; i < 121; i++){
51    i1Sig[i]=0;
52    i2Sig[i]=0;
53    i3Sig[i]=0;
54    i4Sig[i]=0;
55    i5Sig[i]=0;
56    i6Sig[i]=0;
57
58    i1Bkg[i]=0;
59    i2Bkg[i]=0;
60    i3Bkg[i]=0;
61    i4Bkg[i]=0;
62    i5Bkg[i]=0;
63    i6Bkg[i]=0;
64  }
65
66  //Make the set of values for cos1 and cos2
67  for (int i = 0; i < 11; i++){
68    cos1V[i] = i*0.1;
69    cos2V[i] = i*0.1;
70  }
71
72  //Make the set of values for MinLep1
73  for (int i = 0; i < 1; i++){
74    dMinV[i] = i + 47;
75  }
76 
77
78  TProfile *h_sigma = new TProfile("sigma", "", 121, 0, 121, 0, 10);
79  TH2F *h_sigma2D   = new TH2F("sigma2D", "", 11, 0, 1, 11, 0, 1);
80
81  for (int i = 0; i < 2; i++){
82    int nevt = t[i]->GetEntries();
83    for (int iev = 0; iev < nevt; iev++) {
84      t[i]->GetEntry(iev);
85     
86      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
87      double HT   = t[i]->GetLeaf("HT")->GetValue();   
88      if (mt > 300 && HT > 600){
89        iTot[i]++;
90      }   
91    }
92  }
93
94  for (int i = 0; i < 2; i++){
95    int nevt = t[i]->GetEntries();
96
97    int icount = 0;
98    for (int iC1 = 0; iC1 < 11; iC1++){
99      double Cut1 = cos1V[iC1];
100      for (int iC2 = 0; iC2 < 11; iC2++){
101        double Cut2 = cos2V[iC2];
102        int k0 = 11*iC1 + iC2;
103        for (int iC3 = 0; iC3 < 1; iC3++){
104          double Cut3 = dMinV[iC3];
105          icount++;
106           
107          int k = 1*k0 + iC3;
108         
109          cout << "file " << i << endl;
110          cout << "iC1 = " << iC1 << "  iC2 = " << iC2 << "  iC3 = " <<  iC3 << endl;
111          cout << "k = " << k << endl;
112          cout << "icount = " << icount << endl;
113          cout << "Cut1 = " << Cut1 << endl;
114          cout << "Cut2 = " << Cut2 << endl;
115          cout << "Cut3 = " << Cut3 << endl << endl;
116         
117          for (int iev = 0; iev < nevt; iev++) {
118            t[i]->GetEntry(iev);
119           
120            double mt   = t[i]->GetLeaf("mtophad")->GetValue();
121            double HT   = t[i]->GetLeaf("HT")->GetValue();
122            double x    = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
123            double y    = t[i]->GetLeaf("cosWqq")->GetValue();
124            double MinLep1 = t[i]->GetLeaf("MinLep1")->GetValue();
125           
126            double cos1 = fabs(x);
127            double cos2 = fabs(y);
128     
129            //For cat cos + cat MinLep < Cut3
130            if (mt > 300 && HT > 600){
131              if (cos1 < Cut1 && cos2 < Cut2 && HT < 800){
132                if (i==0) i1Sig[k]++;
133                if (i==1) i1Bkg[k]++;
134              }
135       
136              if (cos1 > Cut1 && cos2 > Cut2 && HT < 800){
137                if (i==0) i2Sig[k]++;
138                if (i==1) i2Bkg[k]++;
139              }
140       
141              if ((cos1 < Cut1 && cos2 > Cut2 && HT < 800) || (cos1 > Cut1 && cos2 < Cut2 && HT < 800)){
142                if (i==0) i3Sig[k]++;
143                if (i==1) i3Bkg[k]++;           
144              }
145       
146              if (cos1 < Cut1 && cos2 < Cut2 && HT > 800){
147                if (i==0) i4Sig[k]++;
148                if (i==1) i4Bkg[k]++;
149              }
150       
151              if (cos1 > Cut1 && cos2 > Cut2 && HT > 800){
152                if (i==0) i5Sig[k]++;
153                if (i==1) i5Bkg[k]++;
154              }
155       
156              if ((cos1 < Cut1 && cos2 > Cut2 && HT > 800) || (cos1 > Cut1 && cos2 < Cut2 && HT > 800)){
157                if (i==0) i6Sig[k]++;
158                if (i==1) i6Bkg[k]++;
159              }
160            }
161          }//Read each file
162        }//End of loop on iC3
163      }//End of loop on iC2
164    }//End of loop on iC1
165  }//End of reading 2 files
166
167  //Now the table contains number of events for different categories is filled
168  double N1_bar = iTot[0]*Wei[0];
169  double N0_bar = iTot[1]*Wei[1];
170
171  cout << "N1_bar = " << N1_bar << endl;
172  cout << "N0_bar = " << N0_bar << endl;
173
174  double sigmaMin = 1e30;
175  int iBest = -100;
176  int i1Best = -100;
177  int i2Best = -100;
178  double sigV[121];
179 
180  for (int i = 0; i < 121; i++){
181    double F1Cat1 = (double)i1Sig[i]/iTot[0];
182    double F1Cat2 = (double)i2Sig[i]/iTot[0];
183    double F1Cat3 = (double)i3Sig[i]/iTot[0];
184    double F1Cat4 = (double)i4Sig[i]/iTot[0];
185    double F1Cat5 = (double)i5Sig[i]/iTot[0];
186    double F1Cat6 = (double)i6Sig[i]/iTot[0];
187
188    double F0Cat1 = (double)i1Bkg[i]/iTot[1];
189    double F0Cat2 = (double)i2Bkg[i]/iTot[1];
190    double F0Cat3 = (double)i3Bkg[i]/iTot[1];
191    double F0Cat4 = (double)i4Bkg[i]/iTot[1];
192    double F0Cat5 = (double)i5Bkg[i]/iTot[1];
193    double F0Cat6 = (double)i6Bkg[i]/iTot[1];
194
195    //Store the value of the best cuts
196    if (i==95){
197      fprintf(fRes3,"%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f",F1Cat1,F1Cat2,F1Cat3,F1Cat4,F1Cat5,F1Cat6,F0Cat1,F0Cat2,F0Cat3,F0Cat4,F0Cat5,F0Cat6);
198      fclose(fRes3);
199    }
200
201    double T1 = F1Cat1*F1Cat1/(N1_bar*F1Cat1 + N0_bar*F0Cat1);
202    if (F1Cat1==0 && F0Cat1==0) T1 = 0;
203    double T2 = F1Cat2*F1Cat2/(N1_bar*F1Cat2 + N0_bar*F0Cat2);
204    if (F1Cat2==0 && F0Cat2==0) T2 = 0;
205    double T3 = F1Cat3*F1Cat3/(N1_bar*F1Cat3 + N0_bar*F0Cat3);
206    if (F1Cat3==0 && F0Cat3==0) T3 = 0;
207    double T4 = F1Cat4*F1Cat4/(N1_bar*F1Cat4 + N0_bar*F0Cat4);
208    if (F1Cat4==0 && F0Cat4==0) T4 = 0;
209    double T5 = F1Cat5*F1Cat5/(N1_bar*F1Cat5 + N0_bar*F0Cat5);
210    if (F1Cat5==0 && F0Cat5==0) T5 = 0;
211    double T6 = F1Cat6*F1Cat6/(N1_bar*F1Cat6 + N0_bar*F0Cat6);
212    if (F1Cat6==0 && F0Cat6==0) T6 = 0;
213
214    double sigmaInv = T1 + T2 + T3 + T4 + T5 + T6;
215    double sigma2 = 1/N1_bar/sigmaInv;
216    double sigma = sqrt(sigma2);
217    sigV[i] = sigma;
218
219    cout << "sigmaInv = " << sigmaInv << " sigma = " << sigma << endl;
220
221    h_sigma->Fill(i,sigma);
222
223    if (sigma < sigmaMin){
224      sigmaMin = sigma;
225      iBest = i;
226    }
227  }
228
229  cout << "iTot[0] = " << iTot[0] << endl;
230  cout << "iTot[1] = " << iTot[1] << endl;
231  cout << "sigmaMin = " << sigmaMin << endl;
232  cout << "iBest = " << iBest << endl;
233
234  for (int i = 0; i < 11; i++){
235    for (int j = 0; j < 11; j++){
236      int iC0 = 11*i + j;
237      for (int k = 0; k < 1; k++){
238        int iC = 1*iC0 + k;
239        double Cut1 = cos1V[i];
240        double Cut2 = cos2V[j];
241
242        h_sigma2D->Fill(cos1V[i],cos2V[j],sigV[iC]);
243
244        if (iC==iBest) {
245          i1Best = i;
246          i2Best = j;
247          cout << "i1Best = " << i << endl;
248          cout << "i2Best = " << j << endl;
249          cout << "cos1 = " << cos1V[i] << endl;
250          cout << "cos2 = " << cos2V[j] << endl;
251          cout << "dMin = " << dMinV[k] << endl;
252        }
253
254        if (cos1V[i]==0.6)
255          cout << "Cut1 = 0.6" << endl;
256        if (cos2V[j]==0.6)
257          cout << "Cut2 = 0.6" << endl;
258
259        if (cos1V[i] == 0.6 && cos2V[j] == 0.6){
260          cout << "sigma when Cut1=0.6 and Cut2=0.6 : " << sigV[iC] << endl;
261        }
262        if (i==6 && j==6){
263          cout << "sigma when i=6 and j=6 : " << sigV[iC] << endl;
264        }
265        if (i==i1Best&&j==i2Best)
266          cout << "sigma when i1Best and i2Best : " << sigV[iC] << endl;;
267      }
268    }
269  }
270
271  TCanvas *c1 = new TCanvas("c1", "", 500,400);
272  c1->cd();
273  h_sigma->SetLineWidth(3);
274  h_sigma->Draw();
275
276  TCanvas *c2 = new TCanvas("c2", "", 500,400);
277  c2->cd();
278  h_sigma2D->Draw("CONT4Z");
279
280  c1->Print("Sigma_WithoutMinLep.ps");
281  c2->Print("Sigma2D_WithoutMinLep.ps");
282}
Note: See TracBrowser for help on using the repository browser.