ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/claudioc/ScansICHEP2012/example2011/makeNtuple.C
Revision: 1.1
Committed: Fri Jun 8 16:55:13 2012 UTC (12 years, 10 months ago) by claudioc
Content type: text/plain
Branch: MAIN
CVS Tags: CommonAN1stRelease, HEAD
Error occurred while calculating annotation data.
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 7 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_sr1.txt > SR1_sorted.txt
17 // sort t1tttt_sr2.txt > SR2_sorted.txt
18 // sort t1tttt_sr3.txt > SR3_sorted.txt
19 // sort t1tttt_sr4.txt > SR4_sorted.txt
20 // sort t1tttt_sr5.txt > SR5_sorted.txt
21 // sort t1tttt_sr6.txt > SR6_sorted.txt
22 // sort t1tttt_sr7.txt > SR7_sorted.txt
23 //
24 // Make sure that the file names, which are
25 // hardwired (sorry) are correct.
26 //
27 // Also reads in the gluino xsection txt file.
28 // Assumes that the gluino mass increases monotonically
29 // from one line to the next
30 //
31 //
32 // Usage:
33 // root> .L makeNtuple()
34 // root> makeNtuple()
35 //-------------------------------------------
36
37 //
38 //
39
40 void makeNtuple(){
41
42 //-----------------------------------------------
43 // The name of the file with the output ntuple
44 //-----------------------------------------------
45 char* ntFile = "ntuple.root";
46
47
48 //-----------------------------------------------
49 // Raed in the gluino pair production xsection
50 //-----------------------------------------------
51 float xsec_mass[2000];
52 float xsec[2000];
53 float xsecup[2000];
54 float xsecdwn[2000];
55 int nxsec=0;
56 ifstream fromfile("gluino_xsec_7Tev.txt");
57 if ( !fromfile.good() ) {
58 cout << "gluino xsec file does not exist" << endl ;
59 return;
60 }
61 cout << "Reading gluino xsection" << endl;
62 int xs=0;
63 while (!fromfile.eof()) {
64 TString d1, d2, d3, d4;
65 fromfile >> d1 >> xsec_mass[xs] >> d2 >> xsec[xs] >> d3
66 >> xsecup[xs] >> d4 >> xsecdwn[xs];
67 xs++;
68 nxsec++;
69 }
70 fromfile.close();
71 nxsec = nxsec-1; //last line is junk
72
73 for (int i=0; i<nxsec; i++) {
74 cout<<xsec_mass[i]<<" "<<xsec[i]<<" "<<xsecup[i]<<" "<< xsecdwn[i]<< endl;
75 }
76
77 //--------------------------------------------------------
78 // Check that the gluino mass in the xsection file increases monotonically
79 //---------------------------------------------------------
80 float previous = xsec_mass[0];
81 for (int i=1; i<nxsec; i++) {
82 if (xsec_mass[i] < previous) {
83 cout << "Mass in xsection file does not increase monotonically..Abort"<<endl;
84 return;
85 }
86 previous = xsec_mass[i];
87 }
88
89 //-----------------------------------------------
90 // Define the variables in the txt files
91 //-----------------------------------------------
92 float glmass[7][1000];
93 float lspmass[7][1000];
94 float eff[7][1000];
95 float err[7][1000];
96 float obslim[7][1000];
97 float explim[7][1000];
98 bool usedflag[7][1000]; // a flag to be used later
99 int nentries[7]; // number of entries in the file
100
101 //-----------------------------------------------
102 // The file names, ordered by signal region...
103 //-----------------------------------------------
104 vector<TString> filenames;
105 filenames.push_back("t1tttt_sr1.txt");
106 filenames.push_back("t1tttt_sr2.txt");
107 filenames.push_back("t1tttt_sr3.txt");
108 filenames.push_back("t1tttt_sr4.txt");
109 filenames.push_back("t1tttt_sr5.txt");
110 filenames.push_back("t1tttt_sr6.txt");
111 filenames.push_back("t1tttt_sr7.txt");
112
113 //-----------------------------------------------
114 // Loop over signal regions and load the arrays from the txt files
115 //-----------------------------------------------
116 for (int ireg=0; ireg<7; ireg++) {
117 nentries[ireg] = 0;
118 ifstream fromfile(filenames.at(ireg));
119
120 // if the file does not exist, go to the next one
121 if ( !fromfile.good() ) {
122 cout << "file " << filenames.at(ireg) << " does not exist" << endl ;
123 continue;
124 }
125 // the file exist, read from it
126 cout << "Reading from " << filenames.at(ireg) << endl;
127 int j=0;
128 while (!fromfile.eof()) {
129 fromfile >> glmass[ireg][j]
130 >> lspmass[ireg][j]
131 >> eff[ireg][j]
132 >> err[ireg][j]
133 >> obslim[ireg][j]
134 >> explim[ireg][j];
135 j++;
136 nentries[ireg]++;
137 }
138 fromfile.close();
139 // the last entry is junk
140 nentries[ireg] = nentries[ireg]-1;
141 }
142 //-------------------------------------------
143 // Now check that the line order is OK
144 // And that the number of lines is consistent
145 //--------------------------------------------
146 bool bad=false;
147 int numbLines = 0;
148 for (int ii=0; ii<7; ii++) {
149 for (int jj=ii+1; jj<7; jj++) {
150 if (nentries[ii]*nentries[jj] !=0 ) {
151 numbLines = nentries[ii];
152 if (nentries[ii] != nentries[jj]) {
153 cout << "Mismatched lines: file " << filenames.at(ii);
154 cout << " has " << nentries[ii] << " lines" << endl;
155 cout << " file " << filenames.at(jj);
156 cout << " has " << nentries[jj] << " lines" << endl;
157 bad=true;
158 }
159 for (int line=0; line<nentries[ii]; line++) {
160 if (glmass[ii][line] != glmass[jj][line]) {
161 cout << "Gluino mass on line " << line
162 << " of files " << filenames.at(ii) << " and "
163 << filenames.at(jj) << " do not match" << endl;
164 bad=true;
165 }
166 if (lspmass[ii][line] != lspmass[jj][line]) {
167 cout << "LSP mass on line " << line
168 << " of files " << filenames.at(ii) << " and "
169 << filenames.at(jj) << " do not match" << endl;
170 bad=true;
171 }
172 }
173 }
174 }
175 }
176 if (bad) return;
177 cout << " The file formats match...Now fill ntuple" <<endl;
178
179 //-----------------------------------------------------
180 // Now the content of the files is loaded in arrays
181 // We are going to stick everything in an ntuple
182 //-----------------------------------------------------
183 TFile* babyFile_ = TFile::Open(Form("%s", ntFile), "RECREATE");
184 babyFile_->cd();
185 babyTree_ = new TTree("tree", "A Baby Ntuple");
186 babyTree_->SetDirectory(0);
187
188 // define the variables
189 float glmass_; // gluino mass
190 float lspmass_; // lsp mass
191 float xsec_; // xsec
192 float xsecup_; // xsec, one sigma high
193 float xsecdwn_; // xsec, one sigma down
194 int bestsr_ = 0; // best signal region
195
196 // Here come the efficiencies for the various SR.
197 // The last one is the "best", based on the best expected limit
198 float effsr1_ = -1.;
199 float effsr2_ = -1.;
200 float effsr3_ = -1.;
201 float effsr4_ = -1.;
202 float effsr5_ = -1.;
203 float effsr6_ = -1.;
204 float effsr7_ = -1.;
205 float effsrb_ = -1.; // b = best SR
206
207 // same as above but for the efficiency error
208 float errsr1_ = -1.;
209 float errsr2_ = -1.;
210 float errsr3_ = -1.;
211 float errsr4_ = -1.;
212 float errsr5_ = -1.;
213 float errsr6_ = -1.;
214 float errsr7_ = -1.;
215 float errsrb_ = -1.; // b = best SR
216
217 // same as above but for observed limit
218 float obslimsr1_ = -1.;
219 float obslimsr2_ = -1.;
220 float obslimsr3_ = -1.;
221 float obslimsr4_ = -1.;
222 float obslimsr5_ = -1.;
223 float obslimsr6_ = -1.;
224 float obslimsr7_ = -1.;
225 float obslimsrb_ = -1.; // b = best SR
226
227 // same as above but for expected limit
228 float explimsr1_ = -1.;
229 float explimsr2_ = -1.;
230 float explimsr3_ = -1.;
231 float explimsr4_ = -1.;
232 float explimsr5_ = -1.;
233 float explimsr6_ = -1.;
234 float explimsr7_ = -1.;
235 float explimsrb_ = -1.; // b = best SR
236
237
238 // put the branches in the tree
239 babyTree_->Branch("glmass", &glmass_);
240 babyTree_->Branch("lspmass", &lspmass_);
241 babyTree_->Branch("xsec", &xsec_);
242 babyTree_->Branch("xsecup", &xsecup_);
243 babyTree_->Branch("xsecdwn", &xsecdwn_);
244 babyTree_->Branch("bestsr", &bestsr_);
245
246
247 babyTree_->Branch("effsr1", &effsr1_);
248 babyTree_->Branch("effsr2", &effsr2_);
249 babyTree_->Branch("effsr3", &effsr3_);
250 babyTree_->Branch("effsr4", &effsr4_);
251 babyTree_->Branch("effsr5", &effsr5_);
252 babyTree_->Branch("effsr6", &effsr6_);
253 babyTree_->Branch("effsr7", &effsr7_);
254 babyTree_->Branch("effsrb", &effsrb_); // b = best SR
255
256 babyTree_->Branch("errsr1", &errsr1_);
257 babyTree_->Branch("errsr2", &errsr2_);
258 babyTree_->Branch("errsr3", &errsr3_);
259 babyTree_->Branch("errsr4", &errsr4_);
260 babyTree_->Branch("errsr5", &errsr5_);
261 babyTree_->Branch("errsr6", &errsr6_);
262 babyTree_->Branch("errsr7", &errsr7_);
263 babyTree_->Branch("errsrb", &errsrb_); // b = best SR
264
265 babyTree_->Branch("obslimsr1", &obslimsr1_);
266 babyTree_->Branch("obslimsr2", &obslimsr2_);
267 babyTree_->Branch("obslimsr3", &obslimsr3_);
268 babyTree_->Branch("obslimsr4", &obslimsr4_);
269 babyTree_->Branch("obslimsr5", &obslimsr5_);
270 babyTree_->Branch("obslimsr6", &obslimsr6_);
271 babyTree_->Branch("obslimsr7", &obslimsr7_);
272 babyTree_->Branch("obslimsrb", &obslimsrb_); // b = best SR
273
274 babyTree_->Branch("explimsr1", &explimsr1_);
275 babyTree_->Branch("explimsr2", &explimsr2_);
276 babyTree_->Branch("explimsr3", &explimsr3_);
277 babyTree_->Branch("explimsr4", &explimsr4_);
278 babyTree_->Branch("explimsr5", &explimsr5_);
279 babyTree_->Branch("explimsr6", &explimsr6_);
280 babyTree_->Branch("explimsr7", &explimsr7_);
281 babyTree_->Branch("explimsrb", &explimsrb_); // b = best SR
282
283 // fill the variables.
284 // Note: the best region is the one with the smallest
285 // ratio of expected limit and efficiency
286 cout << "Filling an ntuple with " << numbLines << " entries" << endl;
287 for (int il=0; il<numbLines; il++) {
288 float best = 999999.;
289 int ibest = -1;
290 if (nentries[0] != 0) {
291 glmass_ = glmass[0][il];
292 lspmass_ = lspmass[0][il];
293 effsr1_ = eff[0][il];
294 errsr1_ = eff[0][il];
295 obslimsr1_ = obslim[0][il];
296 explimsr1_ = explim[0][il];
297 float temp = explimsr1_/effsr1_;
298 if (temp < best) {
299 best = temp;
300 ibest = 1;
301 }
302 }
303 if (nentries[1] != 0) {
304 glmass_ = glmass[1][il];
305 lspmass_ = lspmass[1][il];
306 effsr2_ = eff[1][il];
307 errsr2_ = eff[1][il];
308 obslimsr2_ = obslim[1][il];
309 explimsr2_ = explim[1][il];
310 float temp = explimsr2_/effsr2_;
311 if (temp < best) {
312 best = temp;
313 ibest = 2;
314 }
315 }
316 if (nentries[2] != 0) {
317 glmass_ = glmass[2][il];
318 lspmass_ = lspmass[2][il];
319 effsr3_ = eff[2][il];
320 errsr3_ = eff[2][il];
321 obslimsr3_ = obslim[2][il];
322 explimsr3_ = explim[2][il];
323 float temp = explimsr3_/effsr3_;
324 if (temp < best) {
325 best = temp;
326 ibest = 3;
327 }
328 }
329 if (nentries[3] != 0) {
330 glmass_ = glmass[3][il];
331 lspmass_ = lspmass[3][il];
332 effsr4_ = eff[3][il];
333 errsr4_ = eff[3][il];
334 obslimsr4_ = obslim[3][il];
335 explimsr4_ = explim[3][il];
336 float temp = explimsr4_/effsr4_;
337 if (temp < best) {
338 best = temp;
339 ibest = 4;
340 }
341 }
342 if (nentries[4] != 0) {
343 glmass_ = glmass[4][il];
344 lspmass_ = lspmass[4][il];
345 effsr5_ = eff[4][il];
346 errsr5_ = eff[4][il];
347 obslimsr5_ = obslim[4][il];
348 explimsr5_ = explim[4][il];
349 float temp = explimsr5_/effsr5_;
350 if (temp < best) {
351 best = temp;
352 ibest = 5;
353 }
354 }
355 if (nentries[5] != 0) {
356 glmass_ = glmass[5][il];
357 lspmass_ = lspmass[5][il];
358 effsr6_ = eff[5][il];
359 errsr6_ = eff[5][il];
360 obslimsr6_ = obslim[5][il];
361 explimsr6_ = explim[5][il];
362 float temp = explimsr6_/effsr6_;
363 if (temp < best) {
364 best = temp;
365 ibest = 6;
366 }
367 }
368 if (nentries[6] != 0) {
369 glmass_ = glmass[6][il];
370 lspmass_ = lspmass[6][il];
371 effsr7_ = eff[6][il];
372 errsr7_ = eff[6][il];
373 obslimsr7_ = obslim[6][il];
374 explimsr7_ = explim[6][il];
375 float temp = explimsr7_/effsr7_;
376 if (temp < best) {
377 best = temp;
378 ibest = 7;
379 }
380 }
381 bestsr_ = ibest;
382 effsrb_ = eff[ibest-1][il];
383 errsrb_ = eff[ibest-1][il];
384 obslimsrb_ = obslim[ibest-1][il];
385 explimsrb_ = explim[ibest-1][il];
386
387 // Fill the cross-section
388 bool done = false;
389 for (int xs=0; xs<nxsec-1; xs++) {
390 if ( (glmass_ >= xsec_mass[xs] && glmass_ < xsec_mass[xs+1]) ||
391 (glmass_ < xsec_mass[xs] && xs==0) ) {
392
393 float dd1 = glmass_ - xsec_mass[xs];
394 float dd2 = xsec_mass[xs+1] - xsec_mass[xs];
395 xsec_ = xsec[xs] + (xsec[xs+1] - xsec[xs]) *(dd1/dd2);
396 xsecup_ = xsecup[xs] + (xsecup[xs+1] - xsecup[xs]) *(dd1/dd2);
397 xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs+1] - xsecdwn[xs])*(dd1/dd2);
398 done = true;
399 break;
400 }
401 }
402 if (!done) {
403 int xs = nxsec - 1;
404 float dd1 = glmass_ - xsec_mass[xs];
405 float dd2 = xsec_mass[xs] - xsec_mass[xs-1];
406 xsec_ = xsec[xs] + (xsec[xs] - xsec[xs-1]) *(dd1/dd2);
407 xsecup_ = xsecup[xs] + (xsecup[xs] - xsecup[xs-1]) *(dd1/dd2);
408 xsecdwn_ = xsecdwn[xs] + (xsecdwn[xs] - xsecdwn[xs-1])*(dd1/dd2);
409 }
410
411
412 // Now fill the ntuple
413 babyTree_->Fill();
414
415 }
416 // Now close the file
417
418 babyFile_->cd();
419 babyTree_->Write();
420 babyFile_->Close();
421
422
423 }
424