source: huonglan/SPLOT/OptimizeCutWith2DeltaM_NoHT.C

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

re-upload sPlot version
File size: 8.9 KB
Line 
1{
2  gROOT->SetStyle("Plain");
3  TString filename = "";
4  TString type[2] = {"Signal","AllBkg"};
5  TFile *file[2];
6  TTree *t[2];
7  double WeiSig = 0.685*1.43*1000/74873;
8  TFile *fAllBkg = TFile::Open("res_AllBkg.root");
9  TTree *tAllBkg = (TTree*)fAllBkg->Get("T");
10  double WeiBkg = tAllBkg->GetWeight();
11  double Wei[2] = {WeiSig,WeiBkg};
12 
13  for (int i = 0; i < 2; i++){
14    filename = "res_"; filename += type[i]; filename += ".root";
15    file[i] = TFile::Open(filename);
16    t[i] = (TTree*)file[i]->Get("T");
17  }
18 
19  int iTot[2];
20  for (int i = 0; i < 2; i++){
21    iTot[i] = 0;
22  }
23 
24  int  i1S[3025];         int i1B[3025]; 
25  int  i2S[3025];         int i2B[3025]; 
26  int  i3S[3025];         int i3B[3025]; 
27  int  i4S[3025];         int i4B[3025]; 
28  int  i5S[3025];         int i5B[3025]; 
29  int  i6S[3025];         int i6B[3025]; 
30  int  i7S[3025];         int i7B[3025]; 
31  int  i8S[3025];         int i8B[3025]; 
32  int  i9S[3025];         int i9B[3025]; 
33  int i10S[3025];         int i10B[3025];
34  int i11S[3025];         int i11B[3025];
35  int i12S[3025];         int i12B[3025];
36 
37  for (int i = 0; i < 3025; i++){
38    i1S[i]  = 0;       i1B[i]  = 0;
39    i2S[i]  = 0;       i2B[i]  = 0;
40    i3S[i]  = 0;       i3B[i]  = 0;
41    i4S[i]  = 0;       i4B[i]  = 0;
42    i5S[i]  = 0;       i5B[i]  = 0;
43    i6S[i]  = 0;       i6B[i]  = 0;
44    i7S[i]  = 0;       i7B[i]  = 0;
45    i8S[i]  = 0;       i8B[i]  = 0;
46    i9S[i]  = 0;       i9B[i]  = 0;
47    i10S[i] = 0;       i10B[i] = 0;
48    i11S[i] = 0;       i11B[i] = 0;
49    i12S[i] = 0;       i12B[i] = 0;
50  }
51 
52  //#########   OPTIMIZE THE CUTS   ##########
53  double cos1V[5];
54  double cos2V[5];
55  double dMLepV[11];
56  double dMHadV[11];
57  double C1Best = -1e3;
58  double C2Best = -1e3;
59  double C3Best = -1e3;
60  double Rmax   = -1e30;
61  double nSig   = -1e3;
62  double nBkg   = -1e3;
63
64  //Make the set of values for cos1 and cos2
65  for (int i = 0; i < 5; i++){
66    cos1V[i] = i*0.1 + 0.5;
67    cos2V[i] = i*0.1 + 0.5;
68  }
69 
70  //Make the set of values for MinLep1
71  for (int i = 0; i < 11; i++){
72    dMLepV[i] = 2*i + 40;
73    dMHadV[i] = 2*i + 40;
74  }
75 
76  TProfile *h_sigma = new TProfile("sigma", "", 121, 0, 121, 0, 10);
77  TH2F *h_sigma2D   = new TH2F("sigma2D", "", 11, 0, 1, 11, 0, 1);
78 
79  //To calculate the number of HyperNasty events
80  for (int i = 0; i < 2; i++){
81    int nevt = t[i]->GetEntries();
82    for (int iev = 0; iev < nevt; iev++) {
83      t[i]->GetEntry(iev);
84     
85      double mt   = t[i]->GetLeaf("mtophad")->GetValue();
86      double HT   = t[i]->GetLeaf("HT")->GetValue();   
87      if (mt > 300 && HT > 600){
88        iTot[i]++;
89      }   
90    }
91  }
92 
93  double N1_bar = iTot[0]*Wei[0];
94  double N0_bar = iTot[1]*Wei[1];
95 
96  cout << "N1_bar = " << N1_bar << endl;
97  cout << "N0_bar = " << N0_bar << endl;
98 
99  //########  To optimize the cuts now
100  for (int i = 0; i < 2; i++){
101    int nevt = t[i]->GetEntries();
102   
103    int icount = 0;
104    for (int iC1 = 0; iC1 < 5; iC1++){
105      double Cut1 = cos1V[iC1];
106      for (int iC2 = 0; iC2 < 5; iC2++){
107        double Cut2 = cos2V[iC2];
108        int k0 = 5*iC1 + iC2;
109        for (int iC3 = 0; iC3 < 11; iC3++){
110          double Cut3 = dMLepV[iC3];
111          int k1 = 11*k0 + iC3;
112          for (int iC4 = 0; iC4 < 11; iC4++){
113            double Cut4 = dMHadV[iC4];
114            int k = 11*k1 + iC4;
115            icount++;
116           
117            cout << "file  " << i << endl;
118            cout << "iC1 = " << iC1 << "  iC2 = " << iC2 << "  iC3 = " <<  iC3 << " iC4 = " << iC4 << endl;
119            cout << "k      = "  << k      << endl;
120            cout << "icount = "  << icount << endl;
121            cout << "Cut1   = "  << Cut1   << endl;
122            cout << "Cut2   = "  << Cut2   << endl;
123            cout << "Cut3   = "  << Cut3   << endl;
124            cout << "Cut4   = "  << Cut4   << endl << endl;
125           
126            for (int iev = 0; iev < nevt; iev++) {
127              t[i]->GetEntry(iev);
128             
129              double mt    = t[i]->GetLeaf("mtophad")->GetValue();
130              double HT    = t[i]->GetLeaf("HT")->GetValue();
131              double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
132              double y     = t[i]->GetLeaf("cosWqq")->GetValue();
133              double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
134              double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
135           
136              double cos1 = fabs(x);
137              double cos2 = fabs(y);
138             
139              if (mt > 300 && HT > 600) {
140               
141                if (cos1<Cut1 && cos2<Cut2){
142                  if ((dMLep<Cut3) && (dMHad<Cut4)){
143                    if (i==0) i1S[k]++;
144                    if (i==1) i1B[k]++;
145                  }
146                  if ((dMLep<Cut3) && (dMHad>Cut4)){
147                    if (i==0) i2S[k]++;
148                    if (i==1) i2B[k]++;
149                  }
150                  if ((dMLep>Cut3) && (dMHad<Cut4)){
151                    if (i==0) i3S[k]++;
152                    if (i==1) i3B[k]++;
153                  }
154                  if ((dMLep>Cut3) && (dMHad>Cut4)){
155                    if (i==0) i4S[k]++;
156                    if (i==1) i4B[k]++;
157                  }
158                }
159               
160                if (cos1>Cut1 && cos2>Cut2){
161                  if ((dMLep<Cut3) && (dMHad<Cut4)){
162                    if (i==0) i5S[k]++;
163                    if (i==1) i5B[k]++;
164                  }
165                  if ((dMLep<Cut3) && (dMHad>Cut4)){
166                    if (i==0) i6S[k]++;
167                    if (i==1) i6B[k]++;
168                  }
169                  if ((dMLep>Cut3) && (dMHad<Cut4)){
170                    if (i==0) i7S[k]++;
171                    if (i==1) i7B[k]++;
172                  }
173                  if ((dMLep>Cut3) && (dMHad>Cut4)){
174                    if (i==0) i8S[k]++;
175                    if (i==1) i8B[k]++;
176                  }
177                }
178               
179                if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
180                  if ((dMLep<Cut3) && (dMHad<Cut4)){
181                    if (i==0) i9S[k]++;
182                    if (i==1) i9B[k]++;
183                  }
184                  if ((dMLep<Cut3) && (dMHad>Cut4)){
185                    if (i==0) i10S[k]++;
186                    if (i==1) i10B[k]++;
187                  }
188                  if ((dMLep>Cut3) && (dMHad<Cut4)){
189                    if (i==0) i11S[k]++;
190                    if (i==1) i11B[k]++;
191                  }
192                  if ((dMLep>Cut3) && (dMHad>Cut4)){
193                    if (i==0) i12S[k]++;
194                    if (i==1) i12B[k]++;
195                  }
196                }
197              }
198            }//Read each file
199          }//End of loop on iC4
200        }//End of loop on iC3
201      }//End of loop on iC2
202    }//End of loop on iC1
203  }//End of reading 2 files
204 
205  //Now the table contains number of events for different categories is filled
206 
207  double sigmaMin = 1e30;
208  int iBest = -100;
209  int i1Best = -100;
210  int i2Best = -100;
211  int i3Best = -100;
212  int i4Best = -100;
213  double sigV[3025];
214  double F1[12];
215  double F0[12];
216 
217  for (int i = 0; i < 3025; i++){
218    F1[0]   = (double)i1S[i]/iTot[0];
219    F1[1]   = (double)i2S[i]/iTot[0];
220    F1[2]   = (double)i3S[i]/iTot[0];
221    F1[3]   = (double)i4S[i]/iTot[0];
222    F1[4]   = (double)i5S[i]/iTot[0];
223    F1[5]   = (double)i6S[i]/iTot[0];
224    F1[6]   = (double)i7S[i]/iTot[0];
225    F1[7]   = (double)i8S[i]/iTot[0];
226    F1[8]   = (double)i9S[i]/iTot[0];
227    F1[9]   = (double)i10S[i]/iTot[0];
228    F1[10]  = (double)i11S[i]/iTot[0];
229    F1[11]  = (double)i12S[i]/iTot[0];
230   
231    F0[0]   = (double)i1B[i]/iTot[1];
232    F0[1]   = (double)i2B[i]/iTot[1];
233    F0[2]   = (double)i3B[i]/iTot[1];
234    F0[3]   = (double)i4B[i]/iTot[1];
235    F0[4]   = (double)i5B[i]/iTot[1];
236    F0[5]   = (double)i6B[i]/iTot[1];
237    F0[6]   = (double)i7B[i]/iTot[1];
238    F0[7]   = (double)i8B[i]/iTot[1];
239    F0[8]   = (double)i9B[i]/iTot[1];
240    F0[9]   = (double)i10B[i]/iTot[1];
241    F0[10]  = (double)i11B[i]/iTot[1];
242    F0[11]  = (double)i12B[i]/iTot[1];
243   
244    //Store the value of the best cuts
245    //     if (i==73){
246    //       fprintf(fRes1,"%10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f %10.6f",F1Cat1P,F1Cat2P,F1Cat3P,F1Cat1G,F1Cat2G,F1Cat3G,F0Cat1P,F0Cat2P,F0Cat3P,F0Cat1G,F0Cat2G,F0Cat3G);
247    //       fclose(fRes1);
248    //     }
249   
250    double sigmaInv = 0;
251    double term[12];
252    for (int j = 0; j < 12; j++){
253      double term[j] = F1[j]*F1[j]/(N1_bar*F1[j] + N0_bar*F0[j]);
254      if (F1[j]==0 && F0[j]==0)  term[j] = 0;
255      sigmaInv += term[j];
256    }
257   
258    double sigma2 = 1/N1_bar/sigmaInv;
259    double sigma = sqrt(sigma2);
260    sigV[i] = sigma;
261   
262    cout << "sigmaInv = " << sigmaInv << " sigma = " << sigma << endl;
263   
264    h_sigma->Fill(i,sigma);
265   
266    if (sigma < sigmaMin){
267      sigmaMin = sigma;
268      iBest = i;
269    }
270  }
271 
272  cout << "iTot[0] = " << iTot[0] << endl;
273  cout << "iTot[1] = " << iTot[1] << endl;
274  cout << "sigmaMin = " << sigmaMin << endl;
275  cout << "iBest = " << iBest << endl;
276 
277  for (int iC1 = 0; iC1 < 5; iC1++){
278    for (int iC2 = 0; iC2 < 5; iC2++){
279      int k0 = 5*iC1 + iC2;
280      for (int iC3 = 0; iC3 < 11; iC3++){
281        int k1 = 11*k0 + iC3;
282        for (int iC4 = 0; iC4 < 11; iC4++){
283          double Cut1 = cos1V[iC1];
284          double Cut2 = cos2V[iC2];
285          double Cut3 = dMLepV[iC3];
286          double Cut4 = dMHadV[iC4];
287          int k = 11*k1 + iC4;
288         
289          h_sigma2D->Fill(cos1V[iC1],cos2V[iC2],sigV[k]);
290         
291          if (k==iBest) {
292            i1Best = iC1;
293            i2Best = iC2;
294            i3Best = iC3;
295            i4Best = iC4;
296            cout << "i1Best = " << iC1 << endl;
297            cout << "i2Best = " << iC2 << endl;
298            cout << "i3Best = " << iC3 << endl;
299            cout << "i4Best = " << iC4 << endl;
300            cout << "cos1   = " << cos1V[iC1] << endl;
301            cout << "cos2   = " << cos2V[iC2] << endl;
302            cout << "dMLep  = " << dMLepV[iC3] << endl;
303            cout << "dMHad  = " << dMHadV[iC4] << endl;
304          }
305         
306          if (iC1==i1Best && iC2==i2Best && iC3==i3Best && iC4==i4Best)
307            cout << "sigma when i1,2,3,4 Best : " << sigV[k] << endl;;
308        }
309      }
310    }
311  }
312 
313  TCanvas *c1 = new TCanvas("c1", "", 500,400);
314  c1->cd();
315  h_sigma->SetLineWidth(3);
316  h_sigma->Draw();
317 
318  TCanvas *c2 = new TCanvas("c2", "", 500,400);
319  c2->cd();
320  h_sigma2D->Draw("CONT4Z");
321}
Note: See TracBrowser for help on using the repository browser.