ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makeNtuple.C
Revision: 1.1
Committed: Sat Jun 16 02:11:41 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
CVS Tags: CommonAN1stRelease
Log Message:
first

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