ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/code/makeNtuple.C
Revision: 1.4
Committed: Mon Jul 2 04:42:05 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +76 -5 lines
Log Message:
lipstick

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 claudioc 1.4 float explimplus[9][1000];
136     float explimminus[9][1000];
137 claudioc 1.1 bool usedflag[9][1000]; // a flag to be used later
138     int nentries[9]; // number of entries in the file
139    
140     //-----------------------------------------------
141     // The file names, ordered by signal region...
142 claudioc 1.2 // Again, have to sort out the file names
143 claudioc 1.1 //-----------------------------------------------
144 claudioc 1.2 char* fileWhere;
145     if (iflag == 1) fileWhere = "../T1tttt/L2_T1tttt";
146     if (iflag == 2) fileWhere = "../glstop/lsp50/L2_lsp50";
147     if (iflag == 3) fileWhere = "../glstop/lsp250/L2_lsp250";
148     if (iflag == 4) fileWhere = "../glsbottom/chi150/L2_chi150";
149     if (iflag == 5) fileWhere = "../glsbottom/chi300/L2_chi300";
150     if (iflag == 6) fileWhere = "../sbottompair/L2_sbottompair";
151 claudioc 1.1 vector<TString> filenames;
152 claudioc 1.2 for (int kk=0; kk<9; kk++) {
153     filenames.push_back(Form("%s_SR%d_CC.txt", fileWhere,kk));
154     cout << "Will read from file " << filenames[kk] << endl;
155     }
156     //filenames.push_back("../sbottompair/L2_sbottompair_SR0_CC.txt");
157     //filenames.push_back("../sbottompair/L2_sbottompair_SR1_CC.txt");
158     //filenames.push_back("../sbottompair/L2_sbottompair_SR2_CC.txt");
159     //filenames.push_back("../sbottompair/L2_sbottompair_SR3_CC.txt");
160     //filenames.push_back("../sbottompair/L2_sbottompair_SR4_CC.txt");
161     //filenames.push_back("../sbottompair/L2_sbottompair_SR5_CC.txt");
162     //filenames.push_back("../sbottompair/L2_sbottompair_SR6_CC.txt");
163     //filenames.push_back("../sbottompair/L2_sbottompair_SR7_CC.txt");
164     //filenames.push_back("../sbottompair/L2_sbottompair_SR8_CC.txt");
165 claudioc 1.1
166    
167     //-----------------------------------------------
168     // Loop over signal regions and load the arrays from the txt files
169     //-----------------------------------------------
170     for (int ireg=0; ireg<9; ireg++) {
171     nentries[ireg] = 0;
172     ifstream fromfile(filenames.at(ireg));
173    
174     // if the file does not exist, go to the next one
175     if ( !fromfile.good() ) {
176     cout << "file " << filenames.at(ireg) << " does not exist" << endl ;
177     continue;
178     }
179     // the file exist, read from it
180     cout << "Reading from " << filenames.at(ireg) << endl;
181     int j=0;
182     while (!fromfile.eof()) {
183     fromfile >> glmass[ireg][j]
184     >> lspmass[ireg][j]
185     >> eff[ireg][j]
186     >> err[ireg][j];
187     j++;
188     nentries[ireg]++;
189     }
190     fromfile.close();
191 claudioc 1.2 // the last entry is junk
192 claudioc 1.1 nentries[ireg] = nentries[ireg]-1;
193     }
194     //-------------------------------------------
195     // Now check that the line order is OK
196     // And that the number of lines is consistent
197     //--------------------------------------------
198     bool bad=false;
199     int numbLines = 0;
200     for (int ii=0; ii<9; ii++) {
201     for (int jj=ii+1; jj<9; jj++) {
202     if (nentries[ii]*nentries[jj] !=0 ) {
203     numbLines = nentries[ii];
204     if (nentries[ii] != nentries[jj]) {
205     cout << "Mismatched lines: file " << filenames.at(ii);
206     cout << " has " << nentries[ii] << " lines" << endl;
207     cout << " file " << filenames.at(jj);
208     cout << " has " << nentries[jj] << " lines" << endl;
209     bad=true;
210     }
211     for (int line=0; line<nentries[ii]; line++) {
212     if (glmass[ii][line] != glmass[jj][line]) {
213     cout << "Gluino mass on line " << line
214     << " of files " << filenames.at(ii) << " and "
215     << filenames.at(jj) << " do not match" << endl;
216     bad=true;
217     }
218     if (lspmass[ii][line] != lspmass[jj][line]) {
219     cout << "LSP mass on line " << line
220     << " of files " << filenames.at(ii) << " and "
221     << filenames.at(jj) << " do not match" << endl;
222     bad=true;
223     }
224     }
225     }
226     }
227     }
228     if (bad) return;
229     cout << " The file formats match...Now fill ntuple" <<endl;
230    
231     //-------------------------------------
232     //Calculate the total uncertainty
233     //--------------------------------------
234     float lumiErr = 0.045;
235     float trgErr = 0.030;
236     float isoErr = 0.100; // 5% per lepton
237     float fixedErr2 = lumiErr*lumiErr + trgErr*trgErr + isoErr*isoErr;
238     for (int il=0; il<numbLines; il++) {
239     for (int sr=0; sr<9; sr++) {
240     err[sr][il] = sqrt(err[sr][il]*err[sr][il] + fixedErr2);
241     }
242     }
243     //-----------------------------------------------------
244     // And now calculate the limits on the fly
245     //----------------------------------------------------
246     for (int sr=0; sr<9; sr++) {
247     cout << "Calculating limits for SR" << sr << endl;
248     for (int il=0; il<numbLines; il++) {
249 claudioc 1.4 float thisErr = err[sr][il];
250     double blahObs = observedLimitOnNumberOfEvents(sr+1, 100.*thisErr, true );
251     double blahExp = expectedLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
252     double blahExpplus = expectedPlusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
253     double blahExpminus = expectedMinusOneSigmaLimitOnNumberOfEvents(sr+1, 100.*thisErr, false );
254     obslim[sr][il] = (float) blahObs;
255     explim[sr][il] = (float) blahExp;
256     explimplus[sr][il] = (float) blahExpplus;
257     explimminus[sr][il] = (float) blahExpminus;
258 claudioc 1.1 }
259     }
260     //-----------------------------------------------------
261     // Now the content of the files is loaded in arrays
262     // We are going to stick everything in an ntuple
263     //-----------------------------------------------------
264     TFile* babyFile_ = TFile::Open(Form("%s", ntFile), "RECREATE");
265     babyFile_->cd();
266     babyTree_ = new TTree("tree", "A Baby Ntuple");
267     babyTree_->SetDirectory(0);
268    
269     // define the variables
270     float glmass_; // gluino mass
271     float lspmass_; // lsp mass
272     float xsec_; // xsec
273     float xsecup_; // xsec, one sigma high
274     float xsecdwn_; // xsec, one sigma down
275     int bestsr_ = 0; // best signal region
276    
277     // Here come the efficiencies for the various SR.
278     // The last one is the "best", based on the best expected limit
279     float effsr0_ = -1.;
280     float effsr1_ = -1.;
281     float effsr2_ = -1.;
282     float effsr3_ = -1.;
283     float effsr4_ = -1.;
284     float effsr5_ = -1.;
285     float effsr6_ = -1.;
286     float effsr7_ = -1.;
287     float effsr8_ = -1.;
288     float effsrb_ = -1.; // b = best SR
289    
290     // same as above but for the efficiency error
291     float errsr0_ = -1.;
292     float errsr1_ = -1.;
293     float errsr2_ = -1.;
294     float errsr3_ = -1.;
295     float errsr4_ = -1.;
296     float errsr5_ = -1.;
297     float errsr6_ = -1.;
298     float errsr7_ = -1.;
299     float errsr8_ = -1.;
300     float errsrb_ = -1.; // b = best SR
301    
302     // same as above but for observed limit
303     float obslimsr0_ = -1.;
304     float obslimsr1_ = -1.;
305     float obslimsr2_ = -1.;
306     float obslimsr3_ = -1.;
307     float obslimsr4_ = -1.;
308     float obslimsr5_ = -1.;
309     float obslimsr6_ = -1.;
310     float obslimsr7_ = -1.;
311     float obslimsr8_ = -1.;
312     float obslimsrb_ = -1.; // b = best SR
313    
314     // same as above but for expected limit
315     float explimsr0_ = -1.;
316     float explimsr1_ = -1.;
317     float explimsr2_ = -1.;
318     float explimsr3_ = -1.;
319     float explimsr4_ = -1.;
320     float explimsr5_ = -1.;
321     float explimsr6_ = -1.;
322     float explimsr7_ = -1.;
323     float explimsr8_ = -1.;
324     float explimsrb_ = -1.; // b = best SR
325    
326 claudioc 1.4 float explimplussr0_ = -1.;
327     float explimplussr1_ = -1.;
328     float explimplussr2_ = -1.;
329     float explimplussr3_ = -1.;
330     float explimplussr4_ = -1.;
331     float explimplussr5_ = -1.;
332     float explimplussr6_ = -1.;
333     float explimplussr7_ = -1.;
334     float explimplussr8_ = -1.;
335     float explimplussrb_ = -1.; // b = best SR
336    
337     float explimminussr0_ = -1.;
338     float explimminussr1_ = -1.;
339     float explimminussr2_ = -1.;
340     float explimminussr3_ = -1.;
341     float explimminussr4_ = -1.;
342     float explimminussr5_ = -1.;
343     float explimminussr6_ = -1.;
344     float explimminussr7_ = -1.;
345     float explimminussr8_ = -1.;
346     float explimminussrb_ = -1.; // b = best SR
347    
348 claudioc 1.1
349     // put the branches in the tree
350     babyTree_->Branch("glmass", &glmass_);
351     babyTree_->Branch("lspmass", &lspmass_);
352     babyTree_->Branch("xsec", &xsec_);
353     babyTree_->Branch("xsecup", &xsecup_);
354     babyTree_->Branch("xsecdwn", &xsecdwn_);
355     babyTree_->Branch("bestsr", &bestsr_);
356    
357    
358     babyTree_->Branch("effsr0", &effsr0_);
359     babyTree_->Branch("effsr1", &effsr1_);
360     babyTree_->Branch("effsr2", &effsr2_);
361     babyTree_->Branch("effsr3", &effsr3_);
362     babyTree_->Branch("effsr4", &effsr4_);
363     babyTree_->Branch("effsr5", &effsr5_);
364     babyTree_->Branch("effsr6", &effsr6_);
365     babyTree_->Branch("effsr7", &effsr7_);
366     babyTree_->Branch("effsr8", &effsr8_);
367     babyTree_->Branch("effsrb", &effsrb_); // b = best SR
368    
369     babyTree_->Branch("errsr0", &errsr0_);
370     babyTree_->Branch("errsr1", &errsr1_);
371     babyTree_->Branch("errsr2", &errsr2_);
372     babyTree_->Branch("errsr3", &errsr3_);
373     babyTree_->Branch("errsr4", &errsr4_);
374     babyTree_->Branch("errsr5", &errsr5_);
375     babyTree_->Branch("errsr6", &errsr6_);
376     babyTree_->Branch("errsr7", &errsr7_);
377     babyTree_->Branch("errsr8", &errsr8_);
378     babyTree_->Branch("errsrb", &errsrb_); // b = best SR
379    
380     babyTree_->Branch("obslimsr0", &obslimsr0_);
381     babyTree_->Branch("obslimsr1", &obslimsr1_);
382     babyTree_->Branch("obslimsr2", &obslimsr2_);
383     babyTree_->Branch("obslimsr3", &obslimsr3_);
384     babyTree_->Branch("obslimsr4", &obslimsr4_);
385     babyTree_->Branch("obslimsr5", &obslimsr5_);
386     babyTree_->Branch("obslimsr6", &obslimsr6_);
387     babyTree_->Branch("obslimsr7", &obslimsr7_);
388     babyTree_->Branch("obslimsr8", &obslimsr8_);
389     babyTree_->Branch("obslimsrb", &obslimsrb_); // b = best SR
390    
391     babyTree_->Branch("explimsr0", &explimsr0_);
392     babyTree_->Branch("explimsr1", &explimsr1_);
393     babyTree_->Branch("explimsr2", &explimsr2_);
394     babyTree_->Branch("explimsr3", &explimsr3_);
395     babyTree_->Branch("explimsr4", &explimsr4_);
396     babyTree_->Branch("explimsr5", &explimsr5_);
397     babyTree_->Branch("explimsr6", &explimsr6_);
398     babyTree_->Branch("explimsr7", &explimsr7_);
399     babyTree_->Branch("explimsr8", &explimsr8_);
400     babyTree_->Branch("explimsrb", &explimsrb_); // b = best SR
401    
402 claudioc 1.4 babyTree_->Branch("explimplussr0", &explimplussr0_);
403     babyTree_->Branch("explimplussr1", &explimplussr1_);
404     babyTree_->Branch("explimplussr2", &explimplussr2_);
405     babyTree_->Branch("explimplussr3", &explimplussr3_);
406     babyTree_->Branch("explimplussr4", &explimplussr4_);
407     babyTree_->Branch("explimplussr5", &explimplussr5_);
408     babyTree_->Branch("explimplussr6", &explimplussr6_);
409     babyTree_->Branch("explimplussr7", &explimplussr7_);
410     babyTree_->Branch("explimplussr8", &explimplussr8_);
411     babyTree_->Branch("explimplussrb", &explimplussrb_); // b = best SR
412    
413     babyTree_->Branch("explimminussr0", &explimminussr0_);
414     babyTree_->Branch("explimminussr1", &explimminussr1_);
415     babyTree_->Branch("explimminussr2", &explimminussr2_);
416     babyTree_->Branch("explimminussr3", &explimminussr3_);
417     babyTree_->Branch("explimminussr4", &explimminussr4_);
418     babyTree_->Branch("explimminussr5", &explimminussr5_);
419     babyTree_->Branch("explimminussr6", &explimminussr6_);
420     babyTree_->Branch("explimminussr7", &explimminussr7_);
421     babyTree_->Branch("explimminussr8", &explimminussr8_);
422     babyTree_->Branch("explimminussrb", &explimminussrb_); // b = best SR
423    
424    
425 claudioc 1.1 // fill the variables.
426     // Note: the best region is the one with the smallest
427     // ratio of expected limit and efficiency
428     cout << "Filling an ntuple with " << numbLines << " entries" << endl;
429 claudioc 1.2 cout << "SR0 is not allowed to be the best region" << endl;
430 claudioc 1.1 for (int il=0; il<numbLines; il++) {
431     float best = 999999.;
432     int ibest = -1;
433    
434     if (nentries[0] != 0) {
435     glmass_ = glmass[0][il];
436     lspmass_ = lspmass[0][il];
437     effsr0_ = eff[0][il];
438     errsr0_ = eff[0][il];
439     obslimsr0_ = obslim[0][il];
440     explimsr0_ = explim[0][il];
441 claudioc 1.4 explimplussr0_ = explimplus[0][il];
442     explimminussr0_ = explimminus[0][il];
443 claudioc 1.1 float temp = explimsr0_/effsr0_;
444 claudioc 1.2 // if (temp < best) {
445     // best = temp;
446     // ibest = 0;
447     // }
448 claudioc 1.1 }
449     if (nentries[1] != 0) {
450     glmass_ = glmass[1][il];
451     lspmass_ = lspmass[1][il];
452     effsr1_ = eff[1][il];
453     errsr1_ = eff[1][il];
454     obslimsr1_ = obslim[1][il];
455     explimsr1_ = explim[1][il];
456 claudioc 1.4 explimplussr1_ = explimplus[1][il];
457     explimminussr1_ = explimminus[1][il];
458 claudioc 1.1 float temp = explimsr1_/effsr1_;
459     if (temp < best) {
460     best = temp;
461     ibest = 1;
462     }
463     }
464     if (nentries[2] != 0) {
465     glmass_ = glmass[2][il];
466     lspmass_ = lspmass[2][il];
467     effsr2_ = eff[2][il];
468     errsr2_ = eff[2][il];
469     obslimsr2_ = obslim[2][il];
470     explimsr2_ = explim[2][il];
471 claudioc 1.4 explimplussr2_ = explimplus[2][il];
472     explimminussr2_ = explimminus[2][il];
473 claudioc 1.1 float temp = explimsr2_/effsr2_;
474     if (temp < best) {
475     best = temp;
476     ibest = 2;
477     }
478     }
479     if (nentries[3] != 0) {
480     glmass_ = glmass[3][il];
481     lspmass_ = lspmass[3][il];
482     effsr3_ = eff[3][il];
483     errsr3_ = eff[3][il];
484     obslimsr3_ = obslim[3][il];
485     explimsr3_ = explim[3][il];
486 claudioc 1.4 explimplussr3_ = explimplus[3][il];
487     explimminussr3_ = explimminus[3][il];
488 claudioc 1.1 float temp = explimsr3_/effsr3_;
489     if (temp < best) {
490     best = temp;
491     ibest = 3;
492     }
493     }
494     if (nentries[4] != 0) {
495     glmass_ = glmass[4][il];
496     lspmass_ = lspmass[4][il];
497     effsr4_ = eff[4][il];
498     errsr4_ = eff[4][il];
499     obslimsr4_ = obslim[4][il];
500     explimsr4_ = explim[4][il];
501 claudioc 1.4 explimplussr4_ = explimplus[4][il];
502     explimminussr4_ = explimminus[4][il];
503 claudioc 1.1 float temp = explimsr4_/effsr4_;
504     if (temp < best) {
505     best = temp;
506     ibest = 4;
507     }
508     }
509     if (nentries[5] != 0) {
510     glmass_ = glmass[5][il];
511     lspmass_ = lspmass[5][il];
512     effsr5_ = eff[5][il];
513     errsr5_ = eff[5][il];
514     obslimsr5_ = obslim[5][il];
515     explimsr5_ = explim[5][il];
516 claudioc 1.4 explimplussr5_ = explimplus[5][il];
517     explimminussr5_ = explimminus[5][il];
518 claudioc 1.1 float temp = explimsr5_/effsr5_;
519     if (temp < best) {
520     best = temp;
521     ibest = 5;
522     }
523     }
524     if (nentries[6] != 0) {
525     glmass_ = glmass[6][il];
526     lspmass_ = lspmass[6][il];
527     effsr6_ = eff[6][il];
528     errsr6_ = eff[6][il];
529     obslimsr6_ = obslim[6][il];
530     explimsr6_ = explim[6][il];
531 claudioc 1.4 explimplussr6_ = explimplus[6][il];
532     explimminussr6_ = explimminus[6][il];
533 claudioc 1.1 float temp = explimsr6_/effsr6_;
534     if (temp < best) {
535     best = temp;
536     ibest = 6;
537     }
538     }
539     if (nentries[7] != 0) {
540     glmass_ = glmass[7][il];
541     lspmass_ = lspmass[7][il];
542     effsr7_ = eff[7][il];
543     errsr7_ = eff[7][il];
544     obslimsr7_ = obslim[7][il];
545     explimsr7_ = explim[7][il];
546 claudioc 1.4 explimplussr7_ = explimplus[7][il];
547     explimminussr7_ = explimminus[7][il];
548 claudioc 1.1 float temp = explimsr7_/effsr7_;
549     if (temp < best) {
550     best = temp;
551     ibest = 7;
552     }
553     }
554     if (nentries[8] != 0) {
555     glmass_ = glmass[8][il];
556     lspmass_ = lspmass[8][il];
557     effsr8_ = eff[8][il];
558     errsr8_ = eff[8][il];
559     obslimsr8_ = obslim[8][il];
560     explimsr8_ = explim[8][il];
561 claudioc 1.4 explimplussr8_ = explimplus[8][il];
562     explimminussr8_ = explimminus[8][il];
563 claudioc 1.1 float temp = explimsr8_/effsr8_;
564     if (temp < best) {
565     best = temp;
566     ibest = 8;
567     }
568     }
569    
570     bestsr_ = ibest;
571     effsrb_ = eff[ibest][il];
572     errsrb_ = eff[ibest][il];
573     obslimsrb_ = obslim[ibest][il];
574     explimsrb_ = explim[ibest][il];
575 claudioc 1.4 explimplussrb_ = explimplus[ibest][il];
576     explimminussrb_ = explimminus[ibest][il];
577 claudioc 1.1
578     // Fill the cross-section
579     bool done = false;
580     for (int xs=0; xs<nxsec-1; xs++) {
581     if ( (glmass_ >= xsec_mass[xs] && glmass_ < xsec_mass[xs+1]) ||
582     (glmass_ < xsec_mass[xs] && xs==0) ) {
583    
584     float dd1 = glmass_ - xsec_mass[xs];
585     float dd2 = xsec_mass[xs+1] - xsec_mass[xs];
586     xsec_ = xsec[xs] + (xsec[xs+1] - xsec[xs]) *(dd1/dd2);
587     xsecup_ = xsecup[xs] + (xsecup[xs+1] - xsecup[xs]) *(dd1/dd2);
588     xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs+1] - xsecdwn[xs])*(dd1/dd2);
589     done = true;
590     break;
591     }
592     }
593     if (!done) {
594     int xs = nxsec - 1;
595     float dd1 = glmass_ - xsec_mass[xs];
596     float dd2 = xsec_mass[xs] - xsec_mass[xs-1];
597     xsec_ = xsec[xs] + (xsec[xs] - xsec[xs-1]) *(dd1/dd2);
598     xsecup_ = xsecup[xs] + (xsecup[xs] - xsecup[xs-1]) *(dd1/dd2);
599     xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs] - xsecdwn[xs-1])*(dd1/dd2);
600     }
601    
602    
603     // Now fill the ntuple
604     babyTree_->Fill();
605    
606     }
607     // Now close the file
608    
609     babyFile_->cd();
610     babyTree_->Write();
611     babyFile_->Close();
612    
613    
614     }
615