ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makeNtuple.C
Revision: 1.3
Committed: Thu Jun 28 06:08:19 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
Changes since 1.2: +1 -1 lines
Log Message:
fixed limits

File Contents

# Content
1 //-----------------------------------------------
2 // Takes the text files and dumps the contents
3 // in an ntuple.
4 // It's OK if a text file does not exist.
5 // However, the text files need to have the
6 // same order of gluino-mass LSP-mass etc...
7 // Meaning that the lines must correspond
8 // If they do not, the program will barf
9 //
10 // It expects that there are files for 9 signal regions.
11 //
12 // The txt files are supposed to be "sorted".
13 // If they are not, you may want to sort them
14 // ahead of time, eg:
15 //
16 // sort t1tttt_sr0.txt > SR0_sorted.txt
17 // sort t1tttt_sr1.txt > SR1_sorted.txt
18 // sort t1tttt_sr2.txt > SR2_sorted.txt
19 // sort t1tttt_sr3.txt > SR3_sorted.txt
20 // sort t1tttt_sr4.txt > SR4_sorted.txt
21 // sort t1tttt_sr5.txt > SR5_sorted.txt
22 // sort t1tttt_sr6.txt > SR6_sorted.txt
23 // sort t1tttt_sr7.txt > SR7_sorted.txt
24 // sort t1tttt_sr8.txt > SR8_sorted.txt
25 //
26 // Make sure that the file names, which are
27 // hardwired (sorry) are correct.
28 //
29 // Also reads in the gluino xsection txt file.
30 // Assumes that the gluino mass increases monotonically
31 // from one line to the next
32 //
33 //
34 // Usage:
35 // root> .L Functions.cc
36 // root> .L makeNtuple()
37 // root> makeNtuple()
38 //-------------------------------------------
39
40 //
41 //
42
43 void makeNtuple(int iflag=0){
44
45 //------------------------------
46 // Check the flags
47 //------------------------------
48 if (iflag < 1 || iflag > 6) {
49 cout << "Usage:" << endl;
50 cout << "makeNtuple(1) for T1tttt" << endl;
51 cout << "makeNtuple(2) for glstop lsp50" << endl;
52 cout << "makeNtuple(3) for glstop lsp250" << endl;
53 cout << "makeNtuple(4) for glsbottom chi150" << endl;
54 cout << "makeNtuple(5) for glsbottom chi300" << endl;
55 cout << "makeNtuple(6) for glsbottom sbpair" << endl;
56 return;
57 }
58 //---------------------------------------------------
59 // This the model name, used for the root file name
60 //----------------------------------------------------
61 if (iflag == 1) modelName = "T1tttt";
62 if (iflag == 2) modelName = "glstop50";
63 if (iflag == 3) modelName = "glstop250";
64 if (iflag == 4) modelName = "glsb150";
65 if (iflag == 5) modelName = "glsb300";
66 if (iflag == 6) modelName = "sbpair";
67
68 //-----------------------------------------------
69 // The name of the file with the output ntuple
70 //-----------------------------------------------
71 char* ntFile = Form("ntuple_%s.root",modelName);
72 cout << "The output ntuple will be: " << ntFile << endl;
73
74 //-----------------------------------------------
75 // Read in the gluino or sbottom pair production xsection
76 //-----------------------------------------------
77 float xsec_mass[2000];
78 float xsec[2000];
79 float xsecup[2000];
80 float xsecdwn[2000];
81 int nxsec=0;
82 char* xsecFile;
83 if (iflag == 6) {
84 xsecFile = "../xsec/xsec_sbottom.txt";
85 } else {
86 xsecFile = "../xsec/xsec_gluino.txt";
87 }
88 ifstream fromfile(Form("%s",xsecFile));
89 if ( !fromfile.good() ) {
90 cout << "gluino xsec file does not exist" << endl ;
91 return;
92 }
93 cout << "Reading xsection from " << xsecFile << endl;
94 // cout << "ASSUMING THAT THE XSEC IS IN fb...TURN IT INTO pb" << endl;
95 int xs=0;
96 while (!fromfile.eof()) {
97 TString d1, d2, d3, d4;
98 fromfile >> d1 >> xsec_mass[xs] >> d2 >> xsec[xs] >> d3
99 >> xsecup[xs] >> d4 >> xsecdwn[xs];
100 // xsec[xs] = xsec[xs]/1000.;
101 // xsecup[xs] = xsecup[xs]/1000.;
102 // xsecdwn[xs] = xsecdwn[xs]/1000.;
103 xs++;
104 nxsec++;
105 }
106 fromfile.close();
107 nxsec = nxsec-1; //last line is junk
108
109
110 //for (int i=0; i<nxsec; i++) {
111 // cout<<xsec_mass[i]<<" "<<xsec[i]<<" "<<xsecup[i]<<" "<< xsecdwn[i]<< endl;
112 //}
113
114 //--------------------------------------------------------
115 // Check that the gluino mass in the xsection file increases monotonically
116 //---------------------------------------------------------
117 float previous = xsec_mass[0];
118 for (int i=1; i<nxsec; i++) {
119 if (xsec_mass[i] < previous) {
120 cout << "Mass in xsection file does not increase monotonically..Abort"<<endl;
121 return;
122 }
123 previous = xsec_mass[i];
124 }
125
126 //-----------------------------------------------
127 // Define the variables in the txt files
128 //-----------------------------------------------
129 float glmass[9][1000];
130 float lspmass[9][1000];
131 float eff[9][1000];
132 float err[9][1000];
133 float obslim[9][1000];
134 float explim[9][1000];
135 bool usedflag[9][1000]; // a flag to be used later
136 int nentries[9]; // number of entries in the file
137
138 //-----------------------------------------------
139 // The file names, ordered by signal region...
140 // Again, have to sort out the file names
141 //-----------------------------------------------
142 char* fileWhere;
143 if (iflag == 1) fileWhere = "../T1tttt/L2_T1tttt";
144 if (iflag == 2) fileWhere = "../glstop/lsp50/L2_lsp50";
145 if (iflag == 3) fileWhere = "../glstop/lsp250/L2_lsp250";
146 if (iflag == 4) fileWhere = "../glsbottom/chi150/L2_chi150";
147 if (iflag == 5) fileWhere = "../glsbottom/chi300/L2_chi300";
148 if (iflag == 6) fileWhere = "../sbottompair/L2_sbottompair";
149 vector<TString> filenames;
150 for (int kk=0; kk<9; kk++) {
151 filenames.push_back(Form("%s_SR%d_CC.txt", fileWhere,kk));
152 cout << "Will read from file " << filenames[kk] << endl;
153 }
154 //filenames.push_back("../sbottompair/L2_sbottompair_SR0_CC.txt");
155 //filenames.push_back("../sbottompair/L2_sbottompair_SR1_CC.txt");
156 //filenames.push_back("../sbottompair/L2_sbottompair_SR2_CC.txt");
157 //filenames.push_back("../sbottompair/L2_sbottompair_SR3_CC.txt");
158 //filenames.push_back("../sbottompair/L2_sbottompair_SR4_CC.txt");
159 //filenames.push_back("../sbottompair/L2_sbottompair_SR5_CC.txt");
160 //filenames.push_back("../sbottompair/L2_sbottompair_SR6_CC.txt");
161 //filenames.push_back("../sbottompair/L2_sbottompair_SR7_CC.txt");
162 //filenames.push_back("../sbottompair/L2_sbottompair_SR8_CC.txt");
163
164
165 //-----------------------------------------------
166 // Loop over signal regions and load the arrays from the txt files
167 //-----------------------------------------------
168 for (int ireg=0; ireg<9; ireg++) {
169 nentries[ireg] = 0;
170 ifstream fromfile(filenames.at(ireg));
171
172 // if the file does not exist, go to the next one
173 if ( !fromfile.good() ) {
174 cout << "file " << filenames.at(ireg) << " does not exist" << endl ;
175 continue;
176 }
177 // the file exist, read from it
178 cout << "Reading from " << filenames.at(ireg) << endl;
179 int j=0;
180 while (!fromfile.eof()) {
181 fromfile >> glmass[ireg][j]
182 >> lspmass[ireg][j]
183 >> eff[ireg][j]
184 >> err[ireg][j];
185 j++;
186 nentries[ireg]++;
187 }
188 fromfile.close();
189 // the last entry is junk
190 nentries[ireg] = nentries[ireg]-1;
191 }
192 //-------------------------------------------
193 // Now check that the line order is OK
194 // And that the number of lines is consistent
195 //--------------------------------------------
196 bool bad=false;
197 int numbLines = 0;
198 for (int ii=0; ii<9; ii++) {
199 for (int jj=ii+1; jj<9; jj++) {
200 if (nentries[ii]*nentries[jj] !=0 ) {
201 numbLines = nentries[ii];
202 if (nentries[ii] != nentries[jj]) {
203 cout << "Mismatched lines: file " << filenames.at(ii);
204 cout << " has " << nentries[ii] << " lines" << endl;
205 cout << " file " << filenames.at(jj);
206 cout << " has " << nentries[jj] << " lines" << endl;
207 bad=true;
208 }
209 for (int line=0; line<nentries[ii]; line++) {
210 if (glmass[ii][line] != glmass[jj][line]) {
211 cout << "Gluino mass on line " << line
212 << " of files " << filenames.at(ii) << " and "
213 << filenames.at(jj) << " do not match" << endl;
214 bad=true;
215 }
216 if (lspmass[ii][line] != lspmass[jj][line]) {
217 cout << "LSP mass on line " << line
218 << " of files " << filenames.at(ii) << " and "
219 << filenames.at(jj) << " do not match" << endl;
220 bad=true;
221 }
222 }
223 }
224 }
225 }
226 if (bad) return;
227 cout << " The file formats match...Now fill ntuple" <<endl;
228
229 //-------------------------------------
230 //Calculate the total uncertainty
231 //--------------------------------------
232 float lumiErr = 0.045;
233 float trgErr = 0.030;
234 float isoErr = 0.100; // 5% per lepton
235 float fixedErr2 = lumiErr*lumiErr + trgErr*trgErr + isoErr*isoErr;
236 for (int il=0; il<numbLines; il++) {
237 for (int sr=0; sr<9; sr++) {
238 err[sr][il] = sqrt(err[sr][il]*err[sr][il] + fixedErr2);
239 }
240 }
241 //-----------------------------------------------------
242 // And now calculate the limits on the fly
243 //----------------------------------------------------
244 for (int sr=0; sr<9; sr++) {
245 cout << "Calculating limits for SR" << sr << endl;
246 for (int il=0; il<numbLines; il++) {
247 float thisErr = err[sr][il];
248 double blahObs = observedLimitOnNumberOfEvents(sr+1, 100.*thisErr, true );
249 double blahExp = expectedLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
250 obslim[sr][il] = (float) blahObs;
251 explim[sr][il] = (float) blahExp;
252 }
253 }
254 //-----------------------------------------------------
255 // Now the content of the files is loaded in arrays
256 // We are going to stick everything in an ntuple
257 //-----------------------------------------------------
258 TFile* babyFile_ = TFile::Open(Form("%s", ntFile), "RECREATE");
259 babyFile_->cd();
260 babyTree_ = new TTree("tree", "A Baby Ntuple");
261 babyTree_->SetDirectory(0);
262
263 // define the variables
264 float glmass_; // gluino mass
265 float lspmass_; // lsp mass
266 float xsec_; // xsec
267 float xsecup_; // xsec, one sigma high
268 float xsecdwn_; // xsec, one sigma down
269 int bestsr_ = 0; // best signal region
270
271 // Here come the efficiencies for the various SR.
272 // The last one is the "best", based on the best expected limit
273 float effsr0_ = -1.;
274 float effsr1_ = -1.;
275 float effsr2_ = -1.;
276 float effsr3_ = -1.;
277 float effsr4_ = -1.;
278 float effsr5_ = -1.;
279 float effsr6_ = -1.;
280 float effsr7_ = -1.;
281 float effsr8_ = -1.;
282 float effsrb_ = -1.; // b = best SR
283
284 // same as above but for the efficiency error
285 float errsr0_ = -1.;
286 float errsr1_ = -1.;
287 float errsr2_ = -1.;
288 float errsr3_ = -1.;
289 float errsr4_ = -1.;
290 float errsr5_ = -1.;
291 float errsr6_ = -1.;
292 float errsr7_ = -1.;
293 float errsr8_ = -1.;
294 float errsrb_ = -1.; // b = best SR
295
296 // same as above but for observed limit
297 float obslimsr0_ = -1.;
298 float obslimsr1_ = -1.;
299 float obslimsr2_ = -1.;
300 float obslimsr3_ = -1.;
301 float obslimsr4_ = -1.;
302 float obslimsr5_ = -1.;
303 float obslimsr6_ = -1.;
304 float obslimsr7_ = -1.;
305 float obslimsr8_ = -1.;
306 float obslimsrb_ = -1.; // b = best SR
307
308 // same as above but for expected limit
309 float explimsr0_ = -1.;
310 float explimsr1_ = -1.;
311 float explimsr2_ = -1.;
312 float explimsr3_ = -1.;
313 float explimsr4_ = -1.;
314 float explimsr5_ = -1.;
315 float explimsr6_ = -1.;
316 float explimsr7_ = -1.;
317 float explimsr8_ = -1.;
318 float explimsrb_ = -1.; // b = best SR
319
320
321 // put the branches in the tree
322 babyTree_->Branch("glmass", &glmass_);
323 babyTree_->Branch("lspmass", &lspmass_);
324 babyTree_->Branch("xsec", &xsec_);
325 babyTree_->Branch("xsecup", &xsecup_);
326 babyTree_->Branch("xsecdwn", &xsecdwn_);
327 babyTree_->Branch("bestsr", &bestsr_);
328
329
330 babyTree_->Branch("effsr0", &effsr0_);
331 babyTree_->Branch("effsr1", &effsr1_);
332 babyTree_->Branch("effsr2", &effsr2_);
333 babyTree_->Branch("effsr3", &effsr3_);
334 babyTree_->Branch("effsr4", &effsr4_);
335 babyTree_->Branch("effsr5", &effsr5_);
336 babyTree_->Branch("effsr6", &effsr6_);
337 babyTree_->Branch("effsr7", &effsr7_);
338 babyTree_->Branch("effsr8", &effsr8_);
339 babyTree_->Branch("effsrb", &effsrb_); // b = best SR
340
341 babyTree_->Branch("errsr0", &errsr0_);
342 babyTree_->Branch("errsr1", &errsr1_);
343 babyTree_->Branch("errsr2", &errsr2_);
344 babyTree_->Branch("errsr3", &errsr3_);
345 babyTree_->Branch("errsr4", &errsr4_);
346 babyTree_->Branch("errsr5", &errsr5_);
347 babyTree_->Branch("errsr6", &errsr6_);
348 babyTree_->Branch("errsr7", &errsr7_);
349 babyTree_->Branch("errsr8", &errsr8_);
350 babyTree_->Branch("errsrb", &errsrb_); // b = best SR
351
352 babyTree_->Branch("obslimsr0", &obslimsr0_);
353 babyTree_->Branch("obslimsr1", &obslimsr1_);
354 babyTree_->Branch("obslimsr2", &obslimsr2_);
355 babyTree_->Branch("obslimsr3", &obslimsr3_);
356 babyTree_->Branch("obslimsr4", &obslimsr4_);
357 babyTree_->Branch("obslimsr5", &obslimsr5_);
358 babyTree_->Branch("obslimsr6", &obslimsr6_);
359 babyTree_->Branch("obslimsr7", &obslimsr7_);
360 babyTree_->Branch("obslimsr8", &obslimsr8_);
361 babyTree_->Branch("obslimsrb", &obslimsrb_); // b = best SR
362
363 babyTree_->Branch("explimsr0", &explimsr0_);
364 babyTree_->Branch("explimsr1", &explimsr1_);
365 babyTree_->Branch("explimsr2", &explimsr2_);
366 babyTree_->Branch("explimsr3", &explimsr3_);
367 babyTree_->Branch("explimsr4", &explimsr4_);
368 babyTree_->Branch("explimsr5", &explimsr5_);
369 babyTree_->Branch("explimsr6", &explimsr6_);
370 babyTree_->Branch("explimsr7", &explimsr7_);
371 babyTree_->Branch("explimsr8", &explimsr8_);
372 babyTree_->Branch("explimsrb", &explimsrb_); // b = best SR
373
374 // fill the variables.
375 // Note: the best region is the one with the smallest
376 // ratio of expected limit and efficiency
377 cout << "Filling an ntuple with " << numbLines << " entries" << endl;
378 cout << "SR0 is not allowed to be the best region" << endl;
379 for (int il=0; il<numbLines; il++) {
380 float best = 999999.;
381 int ibest = -1;
382
383 if (nentries[0] != 0) {
384 glmass_ = glmass[0][il];
385 lspmass_ = lspmass[0][il];
386 effsr0_ = eff[0][il];
387 errsr0_ = eff[0][il];
388 obslimsr0_ = obslim[0][il];
389 explimsr0_ = explim[0][il];
390 float temp = explimsr0_/effsr0_;
391 // if (temp < best) {
392 // best = temp;
393 // ibest = 0;
394 // }
395 }
396 if (nentries[1] != 0) {
397 glmass_ = glmass[1][il];
398 lspmass_ = lspmass[1][il];
399 effsr1_ = eff[1][il];
400 errsr1_ = eff[1][il];
401 obslimsr1_ = obslim[1][il];
402 explimsr1_ = explim[1][il];
403 float temp = explimsr1_/effsr1_;
404 if (temp < best) {
405 best = temp;
406 ibest = 1;
407 }
408 }
409 if (nentries[2] != 0) {
410 glmass_ = glmass[2][il];
411 lspmass_ = lspmass[2][il];
412 effsr2_ = eff[2][il];
413 errsr2_ = eff[2][il];
414 obslimsr2_ = obslim[2][il];
415 explimsr2_ = explim[2][il];
416 float temp = explimsr2_/effsr2_;
417 if (temp < best) {
418 best = temp;
419 ibest = 2;
420 }
421 }
422 if (nentries[3] != 0) {
423 glmass_ = glmass[3][il];
424 lspmass_ = lspmass[3][il];
425 effsr3_ = eff[3][il];
426 errsr3_ = eff[3][il];
427 obslimsr3_ = obslim[3][il];
428 explimsr3_ = explim[3][il];
429 float temp = explimsr3_/effsr3_;
430 if (temp < best) {
431 best = temp;
432 ibest = 3;
433 }
434 }
435 if (nentries[4] != 0) {
436 glmass_ = glmass[4][il];
437 lspmass_ = lspmass[4][il];
438 effsr4_ = eff[4][il];
439 errsr4_ = eff[4][il];
440 obslimsr4_ = obslim[4][il];
441 explimsr4_ = explim[4][il];
442 float temp = explimsr4_/effsr4_;
443 if (temp < best) {
444 best = temp;
445 ibest = 4;
446 }
447 }
448 if (nentries[5] != 0) {
449 glmass_ = glmass[5][il];
450 lspmass_ = lspmass[5][il];
451 effsr5_ = eff[5][il];
452 errsr5_ = eff[5][il];
453 obslimsr5_ = obslim[5][il];
454 explimsr5_ = explim[5][il];
455 float temp = explimsr5_/effsr5_;
456 if (temp < best) {
457 best = temp;
458 ibest = 5;
459 }
460 }
461 if (nentries[6] != 0) {
462 glmass_ = glmass[6][il];
463 lspmass_ = lspmass[6][il];
464 effsr6_ = eff[6][il];
465 errsr6_ = eff[6][il];
466 obslimsr6_ = obslim[6][il];
467 explimsr6_ = explim[6][il];
468 float temp = explimsr6_/effsr6_;
469 if (temp < best) {
470 best = temp;
471 ibest = 6;
472 }
473 }
474 if (nentries[7] != 0) {
475 glmass_ = glmass[7][il];
476 lspmass_ = lspmass[7][il];
477 effsr7_ = eff[7][il];
478 errsr7_ = eff[7][il];
479 obslimsr7_ = obslim[7][il];
480 explimsr7_ = explim[7][il];
481 float temp = explimsr7_/effsr7_;
482 if (temp < best) {
483 best = temp;
484 ibest = 7;
485 }
486 }
487 if (nentries[8] != 0) {
488 glmass_ = glmass[8][il];
489 lspmass_ = lspmass[8][il];
490 effsr8_ = eff[8][il];
491 errsr8_ = eff[8][il];
492 obslimsr8_ = obslim[8][il];
493 explimsr8_ = explim[8][il];
494 float temp = explimsr8_/effsr8_;
495 if (temp < best) {
496 best = temp;
497 ibest = 8;
498 }
499 }
500
501 bestsr_ = ibest;
502 effsrb_ = eff[ibest][il];
503 errsrb_ = eff[ibest][il];
504 obslimsrb_ = obslim[ibest][il];
505 explimsrb_ = explim[ibest][il];
506
507 // Fill the cross-section
508 bool done = false;
509 for (int xs=0; xs<nxsec-1; xs++) {
510 if ( (glmass_ >= xsec_mass[xs] && glmass_ < xsec_mass[xs+1]) ||
511 (glmass_ < xsec_mass[xs] && xs==0) ) {
512
513 float dd1 = glmass_ - xsec_mass[xs];
514 float dd2 = xsec_mass[xs+1] - xsec_mass[xs];
515 xsec_ = xsec[xs] + (xsec[xs+1] - xsec[xs]) *(dd1/dd2);
516 xsecup_ = xsecup[xs] + (xsecup[xs+1] - xsecup[xs]) *(dd1/dd2);
517 xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs+1] - xsecdwn[xs])*(dd1/dd2);
518 done = true;
519 break;
520 }
521 }
522 if (!done) {
523 int xs = nxsec - 1;
524 float dd1 = glmass_ - xsec_mass[xs];
525 float dd2 = xsec_mass[xs] - xsec_mass[xs-1];
526 xsec_ = xsec[xs] + (xsec[xs] - xsec[xs-1]) *(dd1/dd2);
527 xsecup_ = xsecup[xs] + (xsecup[xs] - xsecup[xs-1]) *(dd1/dd2);
528 xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs] - xsecdwn[xs-1])*(dd1/dd2);
529 }
530
531
532 // Now fill the ntuple
533 babyTree_->Fill();
534
535 }
536 // Now close the file
537
538 babyFile_->cd();
539 babyTree_->Write();
540 babyFile_->Close();
541
542
543 }
544