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

# 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     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