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

# User Rev Content
1 claudioc 1.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 claudioc 1.2 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 claudioc 1.1
68     //-----------------------------------------------
69     // The name of the file with the output ntuple
70     //-----------------------------------------------
71 claudioc 1.2 char* ntFile = Form("ntuple_%s.root",modelName);
72 claudioc 1.1 cout << "The output ntuple will be: " << ntFile << endl;
73    
74     //-----------------------------------------------
75 claudioc 1.2 // Read in the gluino or sbottom pair production xsection
76 claudioc 1.1 //-----------------------------------------------
77     float xsec_mass[2000];
78     float xsec[2000];
79     float xsecup[2000];
80     float xsecdwn[2000];
81     int nxsec=0;
82 claudioc 1.2 char* xsecFile;
83     if (iflag == 6) {
84     xsecFile = "../xsec/xsec_sbottom.txt";
85     } else {
86     xsecFile = "../xsec/xsec_gluino.txt";
87     }
88 claudioc 1.1 ifstream fromfile(Form("%s",xsecFile));
89     if ( !fromfile.good() ) {
90     cout << "gluino xsec file does not exist" << endl ;
91     return;
92     }
93 claudioc 1.2 cout << "Reading xsection from " << xsecFile << endl;
94 claudioc 1.1 // 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 claudioc 1.2 // Again, have to sort out the file names
141 claudioc 1.1 //-----------------------------------------------
142 claudioc 1.2 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 claudioc 1.1 vector<TString> filenames;
150 claudioc 1.2 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 claudioc 1.1
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 claudioc 1.2 // the last entry is junk
190 claudioc 1.1 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 claudioc 1.3 double blahObs = observedLimitOnNumberOfEvents(sr+1, 100.*thisErr, true );
249 claudioc 1.1 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 claudioc 1.2 cout << "SR0 is not allowed to be the best region" << endl;
379 claudioc 1.1 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 claudioc 1.2 // if (temp < best) {
392     // best = temp;
393     // ibest = 0;
394     // }
395 claudioc 1.1 }
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