source: huonglan/SPLOT/OptimizeNastyCutWith2DeltaM.C

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

re-upload sPlot version
File size: 7.8 KB
Line 
1{
2  gROOT->SetStyle("Plain");
3  TString filename = "";
4  TString type[2] = {"AllBkg","Signal"};
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] = {WeiBkg,WeiSig};
12 
13  //FILE *fRes1 = fopen("CatEfficiency_WithMinLep_NoTouchingHT.txt","w+");
14
15  for (int i = 0; i < 2; i++){
16    filename = "res_"; filename += type[i]; filename += ".root";
17    file[i] = TFile::Open(filename);
18    t[i] = (TTree*)file[i]->Get("T");
19  }
20
21  int iTotS[861];        int iTotB[861];
22  int  i1S[861];         int  i1B[861]; 
23  int  i2S[861];         int  i2B[861]; 
24  int  i3S[861];         int  i3B[861]; 
25  int  i4S[861];         int  i4B[861]; 
26  int  i5S[861];         int  i5B[861]; 
27  int  i6S[861];         int  i6B[861]; 
28  int  i7S[861];         int  i7B[861]; 
29  int  i8S[861];         int  i8B[861]; 
30  int  i9S[861];         int  i9B[861]; 
31  int i10S[861];         int i10B[861];
32  int i11S[861];         int i11B[861];
33  int i12S[861];         int i12B[861];
34
35  for (int i = 0; i < 861; i++){
36    iTotS[i] = 0;      iTotB[i] = 0;
37    i1S[i]  = 0;       i1B[i]  = 0;
38    i2S[i]  = 0;       i2B[i]  = 0;
39    i3S[i]  = 0;       i3B[i]  = 0;
40    i4S[i]  = 0;       i4B[i]  = 0;
41    i5S[i]  = 0;       i5B[i]  = 0;
42    i6S[i]  = 0;       i6B[i]  = 0;
43    i7S[i]  = 0;       i7B[i]  = 0;
44    i8S[i]  = 0;       i8B[i]  = 0;
45    i9S[i]  = 0;       i9B[i]  = 0;
46    i10S[i] = 0;       i10B[i] = 0;
47    i11S[i] = 0;       i11B[i] = 0;
48    i12S[i] = 0;       i12B[i] = 0;
49  }
50
51  //#########   OPTIMIZE THE CUTS   ##########
52  //Set the cut values for cos1, cos1, dMLep, dMHad
53  double Cut1 = 0.6;
54  double Cut2 = 0.7;
55  double Cut3 = 60;
56  double Cut4 = 44;
57  double mtV[21];
58  double htV[41];
59
60  //Make the set of values for mjjj and HT
61  for (int i = 0; i < 21; i++){
62    mtV[i] = 200 + i*5;
63  }
64  for (int i = 0; i < 41; i++){
65    htV[i] = 200 + i*10;
66  }
67
68  TProfile *h_sigma = new TProfile("sigma", "", 861,   0, 861,       0,  10);
69  TH2F *h_sigma2D   = new TH2F("sigma2D",   "",  20, 200, 300, 40, 200, 600);
70
71  //########  Optimize the cuts now
72  int icount = 0;
73  for (int i = 0; i < 2; i++){
74    int nevt = t[i]->GetEntries();
75   
76    for (int iC5 = 0; iC5 < 21; iC5++){
77      double Cut5 = mtV[iC5];
78      for (int iC6 = 0; iC6 < 41; iC6++){
79        icount++;
80        double Cut6 = htV[iC6];
81        int k = 41*iC5 + iC6;
82
83        cout << "file  " << type[i] << endl;
84        cout << "iC5 = " << iC5 << " iC6 = " << iC6 << endl;
85        cout << "k      = "  << k      << endl;
86        cout << "icount = "  << icount << endl;
87        cout << "Cut5   = "  << Cut5   << endl;
88        cout << "Cut6   = "  << Cut6   << endl << endl;
89         
90
91        for (int iev = 0; iev < nevt; iev++) {
92          t[i]->GetEntry(iev);
93           
94          double mt    = t[i]->GetLeaf("mtophad")->GetValue();
95          double HT    = t[i]->GetLeaf("HT")->GetValue();
96          double x     = t[i]->GetLeaf("cos_tlepCMtpair")->GetValue();
97          double y     = t[i]->GetLeaf("cosWqq")->GetValue();
98          double dMLep = t[i]->GetLeaf("MinLep1")->GetValue();
99          double dMHad = t[i]->GetLeaf("MinHad1")->GetValue();
100           
101          double cos1 = fabs(x);
102          double cos2 = fabs(y);
103
104          if (mt > Cut5 && HT > Cut6) {
105            if (i==0) iTotB[k]++;
106            if (i==1) iTotS[k]++;
107
108            if (cos1<Cut1 && cos2<Cut2){
109              if ((dMLep<Cut3) && (dMHad<Cut4)){
110                if (i==0) i1B[k]++;
111                if (i==1) i1S[k]++;
112              }
113              if ((dMLep<Cut3) && (dMHad>Cut4)){
114                if (i==0) i2B[k]++;
115                if (i==1) i2S[k]++;
116              }
117              if ((dMLep>Cut3) && (dMHad<Cut4)){
118                if (i==0) i3B[k]++;
119                if (i==1) i3S[k]++;
120              }
121              if ((dMLep>Cut3) && (dMHad>Cut4)){
122                if (i==0) i4B[k]++;
123                if (i==1) i4S[k]++;
124              }
125            }
126           
127            if (cos1>Cut1 && cos2>Cut2){
128              if ((dMLep<Cut3) && (dMHad<Cut4)){
129                if (i==0) i5B[k]++;
130                if (i==1) i5S[k]++;
131              }
132              if ((dMLep<Cut3) && (dMHad>Cut4)){
133                if (i==0) i6B[k]++;
134                if (i==1) i6S[k]++;
135              }
136              if ((dMLep>Cut3) && (dMHad<Cut4)){
137                if (i==0) i7B[k]++;
138                if (i==1) i7S[k]++;
139              }
140              if ((dMLep>Cut3) && (dMHad>Cut4)){
141                if (i==0) i8B[k]++;
142                if (i==1) i8S[k]++;
143              }
144            }
145           
146            if ((cos1>Cut1 && cos2<Cut2) || (cos1<Cut1 && cos2>Cut2)){
147              if ((dMLep<Cut3) && (dMHad<Cut4)){
148                if (i==0) i9B[k]++;
149                if (i==1) i9S[k]++;
150              }
151              if ((dMLep<Cut3) && (dMHad>Cut4)){
152                if (i==0) i10B[k]++;
153                if (i==1) i10S[k]++;
154              }
155              if ((dMLep>Cut3) && (dMHad<Cut4)){
156                if (i==0) i11B[k]++;
157                if (i==1) i11S[k]++;
158              }
159              if ((dMLep>Cut3) && (dMHad>Cut4)){
160                if (i==0) i12B[k]++;
161                if (i==1) i12S[k]++;
162              }
163            }
164          }
165        }//Read each file
166      }//End of loop on iC5
167    }//End of loop on iC6
168  }//End of reading 2 files
169
170  //Now the table contains number of events for different categories is filled
171
172  double sigmaMin = 1e30;
173  int iBest  = -100;
174  int i5Best = -100;
175  int i6Best = -100;
176  double sigV[861];
177
178  for (int i = 0; i < 861; i++){
179    double F1[12];
180    double F0[12];
181    for (int iF = 0; iF < 12; iF++){
182      F1[iF] = 0;
183      F0[iF] = 0;
184    }
185
186    double N0_bar = iTotB[i]*Wei[0];
187    double N1_bar = iTotS[i]*Wei[1];
188    cout << "i = " << i << endl;
189    cout << "Background iTotB = " << iTotB[i] << endl;
190    cout << "Signal     iTotS = " << iTotS[i] << endl;
191    cout << "N1_bar = " << N1_bar << endl;
192    cout << "N0_bar = " << N0_bar << endl;
193
194    F1[0]   = (double)i1S[i]/iTotS[i];
195    F1[1]   = (double)i2S[i]/iTotS[i];
196    F1[2]   = (double)i3S[i]/iTotS[i];
197    F1[3]   = (double)i4S[i]/iTotS[i];
198    F1[4]   = (double)i5S[i]/iTotS[i];
199    F1[5]   = (double)i6S[i]/iTotS[i];
200    F1[6]   = (double)i7S[i]/iTotS[i];
201    F1[7]   = (double)i8S[i]/iTotS[i];
202    F1[8]   = (double)i9S[i]/iTotS[i];
203    F1[9]   = (double)i10S[i]/iTotS[i];
204    F1[10]  = (double)i11S[i]/iTotS[i];
205    F1[11]  = (double)i12S[i]/iTotS[i];
206
207    F0[0]   = (double)i1B[i]/iTotB[i];
208    F0[1]   = (double)i2B[i]/iTotB[i];
209    F0[2]   = (double)i3B[i]/iTotB[i];
210    F0[3]   = (double)i4B[i]/iTotB[i];
211    F0[4]   = (double)i5B[i]/iTotB[i];
212    F0[5]   = (double)i6B[i]/iTotB[i];
213    F0[6]   = (double)i7B[i]/iTotB[i];
214    F0[7]   = (double)i8B[i]/iTotB[i];
215    F0[8]   = (double)i9B[i]/iTotB[i];
216    F0[9]   = (double)i10B[i]/iTotB[i];
217    F0[10]  = (double)i11B[i]/iTotB[i];
218    F0[11]  = (double)i12B[i]/iTotB[i];
219
220    double sumeffS = 0;
221    double sumeffB = 0;
222    for (int iF = 0; iF < 12; iF++){
223      cout << "F0[" << iF << "] = " << F0[iF] << endl;
224      sumeffB += F0[iF];
225    }
226    for (int iF = 0; iF < 12; iF++){
227      cout << "F1[" << iF << "] = " << F1[iF] << endl;
228      sumeffS += F1[iF];
229    }
230    cout << "sumeffS = " << sumeffS << endl;
231    cout << "sumeffB = " << sumeffB << endl;
232
233    double sigmaInv = 0;
234    double term[12];
235    for (int j = 0; j < 12; j++){
236      double term[j] = F1[j]*F1[j]/(N1_bar*F1[j] + N0_bar*F0[j]);
237      if (F1[j]==0 && F0[j]==0)  term[j] = 0;
238      sigmaInv += term[j];
239    }
240
241    double sigma2 = 1/N1_bar/N1_bar/sigmaInv;
242    double sigma = sqrt(sigma2);
243    sigV[i] = sigma;
244
245    cout << "sigmaInv = " << sigmaInv << " sigma = " << sigma << endl;
246
247    h_sigma->Fill(i,sigma);
248
249    if (sigma < sigmaMin){
250      sigmaMin = sigma;
251      iBest = i;
252    }
253  }
254
255  cout << "sigmaMin = " << sigmaMin << endl;
256  cout << "iBest = " << iBest << endl;
257
258  for (int iC5 = 0; iC5 < 21; iC5++){
259    double Cut5 = mtV[iC5];
260    for (int iC6 = 0; iC6 < 41; iC6++){
261      icount++;
262      double Cut6 = htV[iC6];
263      int k = 41*iC5 + iC6;
264         
265      h_sigma2D->Fill(mtV[iC5],htV[iC6],sigV[k]);
266     
267      if (k==iBest) {
268        i5Best = iC5;
269        i6Best = iC6;
270        cout << "i5Best = " << iC5 << endl;
271        cout << "i6Best = " << iC6 << endl;
272        cout << "mjjj = " << mtV[iC5] << endl;
273        cout << "HT   = " << htV[iC6] << endl;
274      }
275     
276      if (iC5==i5Best && iC6==i6Best)
277        cout << "sigma when i1,2,3,4 Best : " << sigV[k] << endl;;
278    }
279  }
280
281  TCanvas *c1 = new TCanvas("c1", "", 500,400);
282  c1->cd();
283  h_sigma->SetLineWidth(3);
284  h_sigma->Draw();
285
286  TCanvas *c2 = new TCanvas("c2", "", 500,400);
287  c2->cd();
288  h_sigma2D->Draw("CONT4Z");
289}
Note: See TracBrowser for help on using the repository browser.