source: huonglan/SPLOT/NEWREALEASE/OptimizeNastyCut_2cosDelMinLep_NewRel.C @ 17

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

upgrade sPlots
File size: 7.7 KB
Line 
1{
2  gROOT->SetStyle("Plain");
3  TString filename = "";
4  //TString type[2] = {"Top_MCAtNLO_PURewei","Signal"};
5  TString type[2] = {"AllBkgTopMCAtNLO","Signal"};//here v1 just contains Top and Wjets
6  //v3 test : contains all bkg, using the variable hfor
7  TFile *file[2];
8  TTree *t[2];
9  TTree *tw[2];
10
11  for (int i = 0; i < 2; i++){
12    filename = "./MERGEFILES/res_"; filename += type[i]; filename += "v3N.root";
13    file[i] = TFile::Open(filename);
14    t[i]    = (TTree*)file[i]->Get("T");
15    tw[i]   = (TTree*)file[i]->Get("Tw");
16    t[i]->AddFriend(tw[i]);
17  }
18
19  double iTotS[3125];   double iTotB[3125];
20  double   i1S[3125];   double   i1B[3125]; 
21  double   i2S[3125];   double   i2B[3125]; 
22  double   i3S[3125];   double   i3B[3125]; 
23  double   i4S[3125];   double   i4B[3125]; 
24  double   i5S[3125];   double   i5B[3125]; 
25  double   i6S[3125];   double   i6B[3125]; 
26
27  for (int i = 0; i < 3125; i++){
28    iTotS[i] = 0;      iTotB[i] = 0;
29    i1S[i]  = 0;       i1B[i]  = 0;
30    i2S[i]  = 0;       i2B[i]  = 0;
31    i3S[i]  = 0;       i3B[i]  = 0;
32    i4S[i]  = 0;       i4B[i]  = 0;
33    i5S[i]  = 0;       i5B[i]  = 0;
34    i6S[i]  = 0;       i6B[i]  = 0;
35  }
36
37  //#########   OPTIMIZE THE CUTS   ##########
38  //Set the cut values for cos1, cos1, dMLep, dMHad
39  double cos1V[5];
40  double cos2V[5];
41  double dMinV[5];
42  double mtV[5];
43  double htV[5];
44
45  for (int i = 0; i < 5; i++){
46    cos1V[i] = i*0.1 + 0.5;
47    cos2V[i] = i*0.1 + 0.5;
48    dMinV[i] = 45 + i;
49    mtV[i]   = 250 + i*10;
50    htV[i]   = 580 + i*10;
51  }
52
53  TProfile *h_sigma = new TProfile("sigma", "", 3125,   0, 3125,       0,  10);
54  TH2F *h_sigma2D   = new TH2F("sigma2D",   "",  20, 200, 300, 40, 200, 600);
55
56  //########  Optimize the cuts now
57  int icount = 0;
58  for (int i = 0; i < 2; i++){
59    int nevt = t[i]->GetEntries();
60   
61    for (int iC1 = 0; iC1 < 5; iC1++){
62      double Cut1 = cos1V[iC1];
63      for (int iC2 = 0; iC2 < 5; iC2++){
64        double Cut2 = cos2V[iC2];
65        int k1 = 5*iC1 + iC2;
66        for (int iC3 = 0; iC3 < 5; iC3++){
67          double Cut3 = dMinV[iC3];
68          int k2 = 5*k1 + iC3;
69          for (int iC4 = 0; iC4 < 5; iC4++){
70            double Cut4 = mtV[iC4];
71            int k3 = 5*k2 + iC4;
72            for (int iC5 = 0; iC5 < 5; iC5++){
73              double Cut5 = htV[iC5];
74              int k = 5*k3 + iC5;
75
76              icount++;
77              cout << "file  " << type[i] << endl;
78              cout << "k      = "  << k      << endl;
79              cout << "icount = "  << icount << endl;
80              cout << "Cut1   = "  << Cut1   << endl;
81              cout << "Cut2   = "  << Cut2   << endl;
82              cout << "Cut3   = "  << Cut3   << endl;
83              cout << "Cut4   = "  << Cut4   << endl;
84              cout << "Cut5   = "  << Cut5   << endl << endl;
85
86              for (int iev = 0; iev < nevt; iev++) {
87                t[i]->GetEntry(iev);
88               
89                double mt     = t[i]->GetLeaf("mtophad")->GetValue();
90                double HT     = t[i]->GetLeaf("HT")->GetValue();
91                double x      = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
92                double y      = t[i]->GetLeaf("cosWqq")->GetValue();
93                double dMLep  = t[i]->GetLeaf("MinLep1")->GetValue();
94                double dMHad  = t[i]->GetLeaf("MinHad1")->GetValue();
95                double WeiEvt = t[i]->GetLeaf("WeiEvent")->GetValue();
96                double Wei1fb = t[i]->GetLeaf("wei")->GetValue();
97                double WeiTot = WeiEvt*Wei1fb;
98               
99                int hasBJ  = t[i]->GetLeaf("hasBadJet")->GetValue();
100                double hf = t[i]->GetLeaf("hfor")->GetValue();
101               
102                if (hf!=4 && hasBJ==0){
103                 
104                  double cos1 = fabs(x);
105                  double cos2 = fabs(y);
106               
107                  if (mt > Cut4 && HT > Cut5) {
108                    if (i==0) iTotB[k] += WeiTot;
109                    if (i==1) iTotS[k] += WeiTot;
110                   
111                    if (cos1<Cut1 && cos2<Cut2){
112                      if ((dMLep<Cut3)){
113                        if (i==0) i1B[k] += WeiTot;
114                        if (i==1) i1S[k] += WeiTot;
115                      }
116                      if ((dMLep>Cut3)){
117                        if (i==0) i2B[k] += WeiTot;
118                        if (i==1) i2S[k] += WeiTot;
119                      }
120                    }
121                   
122                    if (cos1>Cut1 && cos2>Cut2){
123                      if ((dMLep<Cut3)){
124                        if (i==0) i3B[k] += WeiTot;
125                        if (i==1) i3S[k] += WeiTot;
126                      }
127                      if ((dMLep>Cut3)){
128                        if (i==0) i4B[k] += WeiTot;
129                        if (i==1) i4S[k] += WeiTot;
130                      }
131                    }
132                   
133                    if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
134                      if ((dMLep<Cut3)){
135                        if (i==0) i5B[k] += WeiTot;
136                        if (i==1) i5S[k] += WeiTot;
137                      }
138                      if ((dMLep>Cut3)){
139                        if (i==0) i6B[k] += WeiTot;
140                        if (i==1) i6S[k] += WeiTot;
141                      }
142                    }
143                  }
144                }//Cut on Nasty
145              }//Read each file
146            }//End of loop on iC5
147          }//End of loop on iC4
148        }//End of loop on iC3
149      }//End of loop on iC2
150    }//End of loop on iC1
151  }//End of reading 2 files
152
153  //Now the table contains number of events for different categories is filled
154
155  double sigmaMin = 1e30;
156  int iBest  = -100;
157  int i1Best = -100;
158  int i2Best = -100;
159  int i3Best = -100;
160  int i4Best = -100;
161  int i5Best = -100;
162  double sigV[3125];
163
164  for (int i = 0; i < 3125; i++){
165    double F1[6];
166    double F0[6];
167    for (int iF = 0; iF < 6; iF++){
168      F1[iF] = 0;
169      F0[iF] = 0;
170    }
171
172    double N0_bar = iTotB[i];
173    double N1_bar = iTotS[i];
174    cout << "i = " << i << endl;
175    cout << "Background iTotB = " << iTotB[i] << endl;
176    cout << "Signal     iTotS = " << iTotS[i] << endl;
177    cout << "N1_bar = " << N1_bar << endl;
178    cout << "N0_bar = " << N0_bar << endl;
179
180    F1[0]   = (double)i1S[i]/iTotS[i];
181    F1[1]   = (double)i2S[i]/iTotS[i];
182    F1[2]   = (double)i3S[i]/iTotS[i];
183    F1[3]   = (double)i4S[i]/iTotS[i];
184    F1[4]   = (double)i5S[i]/iTotS[i];
185    F1[5]   = (double)i6S[i]/iTotS[i];
186
187    F0[0]   = (double)i1B[i]/iTotB[i];
188    F0[1]   = (double)i2B[i]/iTotB[i];
189    F0[2]   = (double)i3B[i]/iTotB[i];
190    F0[3]   = (double)i4B[i]/iTotB[i];
191    F0[4]   = (double)i5B[i]/iTotB[i];
192    F0[5]   = (double)i6B[i]/iTotB[i];
193
194    double sumeffS = 0;
195    double sumeffB = 0;
196    for (int iF = 0; iF < 6; iF++){
197      cout << "F0[" << iF << "] = " << F0[iF] << endl;
198      sumeffB += F0[iF];
199    }
200    for (int iF = 0; iF < 6; iF++){
201      cout << "F1[" << iF << "] = " << F1[iF] << endl;
202      sumeffS += F1[iF];
203    }
204    cout << "sumeffS = " << sumeffS << endl;
205    cout << "sumeffB = " << sumeffB << endl;
206
207    double sigmaInv = 0;
208    double term[6];
209    for (int j = 0; j < 6; j++){
210      double term[j] = F1[j]*F1[j]/(N1_bar*F1[j] + N0_bar*F0[j]);
211      if (F1[j]==0 && F0[j]==0)  term[j] = 0;
212      sigmaInv += term[j];
213    }
214
215    double sigma2 = 1/N1_bar/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 << "sigmaMin = " << sigmaMin << endl;
230  cout << "iBest = " << iBest << endl;
231
232  for (int iC1 = 0; iC1 < 5; iC1++){
233    double Cut1 = cos1V[iC1];
234    for (int iC2 = 0; iC2 < 5; iC2++){
235      double Cut2 = cos2V[iC2];
236      int k1 = 5*iC1 + iC2;
237      for (int iC3 = 0; iC3 < 5; iC3++){
238        double Cut3 = dMinV[iC3];
239        int k2 = 5*k1 + iC3;
240        for (int iC4 = 0; iC4 < 5; iC4++){
241          double Cut4 = mtV[iC4];
242          int k3 = 5*k2 + iC4;
243          for (int iC5 = 0; iC5 < 5; iC5++){
244            double Cut5 = htV[iC5];
245            int k = 5*k3 + iC5;
246     
247            if (k==iBest) {
248              i1Best = iC1;
249              i2Best = iC2;
250              i3Best = iC3;
251              i4Best = iC4;
252              i5Best = iC5;
253              cout << "i1Best = " << iC1 << endl;
254              cout << "i2Best = " << iC2 << endl;
255              cout << "i3Best = " << iC3 << endl;
256              cout << "i4Best = " << iC4 << endl;
257              cout << "i5Best = " << iC5 << endl;
258              cout << "cos1Best = " << cos1V[iC1] << endl;
259              cout << "cos2Best = " << cos2V[iC2] << endl;
260              cout << "dMinBest = " << dMinV[iC3] << endl;
261              cout << "mtBest   = " << mtV[iC4] << endl;
262              cout << "htBest   = " << htV[iC5] << endl;
263            }
264          }
265        }
266      }
267    }
268  }
269
270  TCanvas *c1 = new TCanvas("c1", "", 500,400);
271  c1->cd();
272  h_sigma->SetLineWidth(3);
273  h_sigma->Draw();
274
275  TCanvas *c2 = new TCanvas("c2", "", 500,400);
276  c2->cd();
277  h_sigma2D->Draw("CONT4Z");
278}
Note: See TracBrowser for help on using the repository browser.