ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/scripts/plotQCDEstimate.py
Revision: 1.1.1.1 (vendor branch)
Committed: Thu Dec 1 16:28:48 2011 UTC (13 years, 5 months ago) by dhidas
Content type: text/x-python
Branch: dhidas, MAIN
CVS Tags: START, HEAD
Changes since 1.1: +0 -0 lines
Log Message:
osu copy modified

File Contents

# Content
1 #
2 # 31 Mar 2010
3 # Plots (signed) Delta = (est-true)/true, for QCD estimate
4 #
5 #---------------------------------------------------------------
6 from __future__ import division
7 from ROOT import *
8
9 n = 4;
10 nrange = 9;
11 plot_average_est = True;
12
13 # free in 12j, fix in 34j
14 constrain_gaus_mean_12j = False;
15 show_free_fit_res_34j = True;
16
17 # landau free
18 #bool show_free_fit_res_34j = true;
19 #const bool constrain_gaus_mean_12j = true;
20
21
22 # gaus mean-constrained in 12j, fix in 34j
23 def plot_qcd_estimate_gaus_mean12j():
24 constrain_gaus_mean_12j = true;
25 show_free_fits_res_34j = false;
26 plot_qcd_estimate( "gaus" );
27
28 def setStyle():
29
30 gROOT.SetStyle( "Plain" );
31 gStyle.SetTextFont( 42 );
32 gStyle.SetLabelFont( 42, "xy" );
33 gStyle.SetFrameFillColor( 0 );
34 gStyle.SetTitleBorderSize( 0 );
35
36 gStyle.SetTitleH( 0.06 );
37 gStyle.SetPadTopMargin( 0.15 );
38
39 gStyle.SetOptFit( 1 );
40 gStyle.SetPalette( 1 );
41
42 gStyle.SetPadTickX( 1 );
43 gStyle.SetPadTickY( 1 );
44
45 def plot_qcd_estimate( func = "gaus" ):
46
47 setStyle();
48
49 if constrain_gaus_mean_12j:
50 show_free_fit_res_34j = False;
51
52 x = [1, 2, 3, 4]
53 xx = [3, 4]
54 x12j = [1, 2]
55 # Gaus (Free-fits)
56
57 result = open( "est_free_%s.txt" % func );
58 #result.open(Form("gaus_mean_0.3to0.6_12j__rb5/est_free_%s.txt",func));
59
60 m = nrange * ( 4 + 2 );
61 est = []
62 estFix = []
63 max_dev = 0;
64 max_dev_12j = 0;
65 max_dev_34j = 0;
66 max_dev_3j = 0;
67 max_dev_4j = 0;
68 max_dev_34j_fix = 0;
69 max_dev_3j_fix = 0;
70 max_dev_4j_fix = 0;
71
72 # Free fits
73 for i in range( 0, 4 * nrange ):
74 est.append( result.readline() )
75 if fabs( est[-1] ) > max_dev:
76 max_dev = fabs( est[-1] );
77
78 if i % 4 <= 1: #1,2j
79 if fabs( est[-1] ) > max_dev_12j:
80 max_dev_12j = fabs( est[-1] );
81 elif i % 4 == 2 : #3j
82 if fabs( est[-1] ) > max_dev_3j :
83 max_dev_3j = fabs( est[-1] );
84 elif i % 4 == 3: #4mj
85 if fabs( est[-1] ) > max_dev_4j :
86 max_dev_4j = fabs( est[-1] );
87
88 if i % 4 == 0:
89 print #first of 4
90 print i % 4 + 1, "j est[", i, "] = ", est[i], " ";
91 if i % 4 <= 1:
92 print "max_dev_12j = ", max_dev_12j; #12j
93 print endl;
94 result.close();
95
96 # Gaus (Constrained-fits in 3,4j)
97 result.open( "est_fix_%s.txt" % func )
98 #result.open(Form("gaus_mean_0.3to0.6_12j__rb5/est_fix_%s.txt",func));
99
100 estFix2 = [];
101 maxDevFix = 0;
102
103 for i in range( 0, 4 * nrange ):
104
105 # read in 3,4j only
106 #print "amtb: i=", i
107 if not constrain_gaus_mean_12j and i % 4 < 2:
108 continue;
109
110 estFix.append( result.readline() )
111 estFix2.append( estFix[-1] );
112
113
114 print i % 4 + 1, "j estFix[", i, "] = " , estFix[-1], " ", endl;
115
116 #print "amtb: i=", i
117
118
119 if fabs( estFix[-1] ) > maxDevFix :
120 maxDevFix = fabs( estFix[-1] );
121 if fabs( estFix[-1] ) > max_dev_34j_fix :
122 max_dev_34j_fix = fabs( estFix[-1] );
123 if i % 4 == 2:#3j
124 if fabs( estFix[-1] ) > max_dev_3j_fix:
125 max_dev_3j_fix = fabs( estFix[-1] );
126 elif i % 4 == 3: #>=4j
127 if fabs( estFix[-1] ) > max_dev_4j_fix:
128 max_dev_4j_fix = fabs( estFix[-1] );
129
130 for k in range( 0, len( estFix2 ) ):
131 print "amtb: estFix2[", k, "]: ", estFix2[k]
132
133 result.close();
134 max_dev_34j_fix = TMath.Max( max_dev_3j_fix, max_dev_4j_fix );
135 max_dev_34j = TMath.Max( max_dev_3j, max_dev_4j );
136
137 print "maxDevFix: ", maxDevFix
138
139
140 print "\nFor all ranges, |max| deviation is ", max_dev
141 print "For 1,2j, |max| deviation is ", max_dev_12j
142 print "For free fit (3j) |max| deviation is ", max_dev_3j
143 print "For free fit (>=4j) |max| deviation is ", max_dev_4j
144 print "For free fit (3,>=4j), |max| deviation is ", max_dev_34j
145 print "For constrained fit (3j), |max| deviation is ", max_dev_3j_fix
146 print "For constrained fit (>=4j), |max| deviation is ", max_dev_4j_fix
147 print "For constrained fit (3,>=4j), |max| deviation is ", max_dev_34j_fix
148
149
150
151 c2 = TCanvas( "c2", "QCD estimates", 600, 600 );
152
153 y[nrange][4];
154
155 index = 0;
156
157 for i in range( 0, nrange ):
158 for j in range( 0, 4 ):
159 y[i][j] = est[index]; #read in
160 index += 1
161 #print "index="<<index<< endl;
162 #print y[i][j]<<endl;;
163 yy12j = [[0, 0] for x in xrange( nrange )]
164 yy34j = [[0, 0] for x in xrange( nrange )]
165 # ix=0;
166 jj = 0;
167 for i in range( 0, nrange ):
168 jj += 1
169 yy12j[i][0] = estFix[jj];#1j
170 jj += 1
171 yy12j[i][1] = estFix[jj];#2j
172 jj += 1
173 yy34j[i][0] = estFix[jj];#3j
174 jj += 1
175 yy34j[i][1] = estFix[jj];#4mj
176
177 #print "index="<<index<< endl;
178
179
180
181
182
183
184 gStyle.SetMarkerSize( 1.7 );
185 gStyle.SetMarkerStyle( 20 );
186 c2.SetTopMargin( 0.1 );
187 c2.SetLeftMargin( 0.12 );
188 c2.SetRightMargin( 0.35 );
189
190 gr1;
191 gr2;
192 gr3;
193 gr4;
194 gr5;
195 gr6;
196 gr7;
197 gr8;
198 gr9;
199 if constrain_gaus_mean_12j == False:
200 if show_free_fit_res_34j:
201 #draw free-fit results for 1-4j
202 gr1 = TGraph( n, x, y[1 - 1] );
203 gr2 = TGraph( n, x, y[2 - 1] );
204 gr3 = TGraph( n, x, y[3 - 1] );
205 gr4 = TGraph( n, x, y[4 - 1] );
206 gr5 = TGraph( n, x, y[5 - 1] );
207 gr6 = TGraph( n, x, y[6 - 1] );
208 gr7 = TGraph( n, x, y[7 - 1] );
209 gr8 = TGraph( n, x, y[8 - 1] );
210 gr9 = TGraph( n, x, y[9 - 1] );
211 else:
212 #draw: free-fit results for 1-2j, OR
213 # gaus-mean-limited in 1,2j
214 gr1 = TGraph( 2, x12j, y[1 - 1] );
215 gr2 = TGraph( 2, x12j, y[2 - 1] );
216 gr3 = TGraph( 2, x12j, y[3 - 1] );
217 gr4 = TGraph( 2, x12j, y[4 - 1] );
218 gr5 = TGraph( 2, x12j, y[5 - 1] );
219 gr6 = TGraph( 2, x12j, y[6 - 1] );
220 gr7 = TGraph( 2, x12j, y[7 - 1] );
221 gr8 = TGraph( 2, x12j, y[8 - 1] );
222 gr9 = TGraph( 2, x12j, y[9 - 1] );
223
224
225 else:
226 gr1 = TGraph( 2, x12j, yy12j[1 - 1] );
227 gr2 = TGraph( 2, x12j, yy12j[2 - 1] );
228 gr3 = TGraph( 2, x12j, yy12j[3 - 1] );
229 gr4 = TGraph( 2, x12j, yy12j[4 - 1] );
230 gr5 = TGraph( 2, x12j, yy12j[5 - 1] );
231 gr6 = TGraph( 2, x12j, yy12j[6 - 1] );
232 gr7 = TGraph( 2, x12j, yy12j[7 - 1] );
233 gr8 = TGraph( 2, x12j, yy12j[8 - 1] );
234 gr9 = TGraph( 2, x12j, yy12j[9 - 1] );
235
236 # constrained
237 gr_1 = TGraph( 2, xx, yy34j[1 - 1] );
238 gr_2 = TGraph( 2, xx, yy34j[2 - 1] );
239 gr_3 = TGraph( 2, xx, yy34j[3 - 1] );
240 gr_4 = TGraph( 2, xx, yy34j[4 - 1] );
241 gr_5 = TGraph( 2, xx, yy34j[5 - 1] );
242 gr_6 = TGraph( 2, xx, yy34j[6 - 1] );
243 gr_7 = TGraph( 2, xx, yy34j[7 - 1] );
244 gr_8 = TGraph( 2, xx, yy34j[8 - 1] );
245 gr_9 = TGraph( 2, xx, yy34j[9 - 1] );
246
247
248 gr1.SetMarkerColor( kGreen + 1 );
249 gr2.SetMarkerColor( kGreen + 2 );
250 gr3.SetMarkerColor( kGreen + 3 );
251 gr4.SetMarkerColor( kAzure + 7 );
252 gr5.SetMarkerColor( kAzure - 3 );
253 gr6.SetMarkerColor( kBlue );
254 gr7.SetMarkerColor( kOrange );
255 gr8.SetMarkerColor( kOrange - 1 );
256 gr9.SetMarkerColor( kOrange - 6 );
257
258 gr_1.SetMarkerColor( kGreen + 1 );
259 gr_2.SetMarkerColor( kGreen + 2 );
260 gr_3.SetMarkerColor( kGreen + 3 );
261 gr_4.SetMarkerColor( kAzure + 7 );
262 gr_5.SetMarkerColor( kAzure - 3 );
263 gr_6.SetMarkerColor( kBlue );
264 gr_7.SetMarkerColor( kOrange );
265 gr_8.SetMarkerColor( kOrange - 1 );
266 gr_9.SetMarkerColor( kOrange - 6 );
267
268
269 gr_1.SetMarkerStyle( 22 );
270 gr_2.SetMarkerStyle( 22 );
271 gr_3.SetMarkerStyle( 22 );
272 gr_4.SetMarkerStyle( 22 );
273 gr_5.SetMarkerStyle( 22 );
274 gr_6.SetMarkerStyle( 22 );
275 gr_7.SetMarkerStyle( 22 );
276 gr_8.SetMarkerStyle( 22 );
277 gr_9.SetMarkerStyle( 22 );
278
279
280
281 # To get desired x range, draw blank histo
282 gStyle.SetTitleW( 0.9 );
283 gStyle.SetTitleH( 0.05 );#?
284 if func == "gaus":
285 h = TH1D( "h", "Variation of QCD estimates with fit range (Gaussian)", 4, 0.5, 4.5 );
286 elif func == "pol3":
287 h = TH1D( "h", "Variation of QCD estimates with fit range (Pol3)", 4, 0.5, 4.5 );
288 elif func == "landau":
289 h = TH1D( "h", "Variation of QCD estimates with fit range (Landau)", 4, 0.5, 4.5 );
290
291
292 h.SetStats( kFALSE ); # no statistics
293 h.Draw();
294 h.SetYTitle( "Deviation = (Est-True)/True" );
295
296 show_range = int( max_dev ) + 1;
297 if constrain_gaus_mean_12j:
298 show_range = int( max_dev ) + 1;
299
300 print "show_range: ", show_range
301
302 h.GetYaxis().SetRangeUser( 0 - show_range, show_range );
303 h.GetXaxis().SetRangeUser( 0.5, 5.5 );
304 h.GetXaxis().SetBinLabel( 1., "1j" );
305 h.GetXaxis().SetBinLabel( 2., "2j" );
306 h.GetXaxis().SetBinLabel( 3., "3j" );
307 h.GetXaxis().SetBinLabel( 4., "#geq4j" );
308 h.GetXaxis().SetLabelSize( 0.07 );
309 h.GetYaxis().SetTitleOffset( 1.3 );
310
311
312 # Free-fits
313 gr1.Draw( "P" );
314 gr2.Draw( "P" ); #to superimpose graphs, do not re-draw axis
315 gr3.Draw( "P" );
316 gr4.Draw( "P" );
317 gr5.Draw( "P" );
318 gr6.Draw( "P" );
319 gr7.Draw( "P" );
320 gr8.Draw( "P" );
321 gr9.Draw( "P" );
322
323 if not show_free_fit_res_34j:
324 # Constrained fits
325 gr_1.Draw( "P" );
326 gr_2.Draw( "P" );
327 gr_3.Draw( "P" );
328 gr_4.Draw( "P" );
329 gr_5.Draw( "P" );
330 gr_6.Draw( "P" );
331 gr_7.Draw( "P" );
332 gr_8.Draw( "P" );
333 gr_9.Draw( "P" );
334
335 c2.SetGrid( 1, 1 );
336
337 leg = TLegend( 0.65, 0.1, 0.98, 0.9 );
338 leg.SetFillColor( 0 );
339 # leg.SetTextFont(62);
340 if func != "gaus" or not constrain_gaus_mean_12j:
341 leg.AddEntry( gr1, "Free: 0.1-1.0", "p" );
342 leg.AddEntry( gr2, "Free: 0.1-1.2", "p" );
343 leg.AddEntry( gr3, "Free: 0.1-1.4", "p" );
344 leg.AddEntry( gr4, "Free: 0.2-1.1", "p" );
345 leg.AddEntry( gr5, "Free: 0.2-1.3", "p" );
346 leg.AddEntry( gr6, "Free: 0.2-1.5", "p" );
347 leg.AddEntry( gr7, "Free: 0.3-1.2", "p" );
348 leg.AddEntry( gr8, "Free: 0.3-1.4", "p" );
349 leg.AddEntry( gr9, "Free: 0.3-1.6", "p" );
350
351 if func == "gaus" and constrain_gaus_mean_12j:
352 leg.AddEntry( gr1, "#mu-constr.: 0.1-1.0", "p" );
353 leg.AddEntry( gr2, "#mu-constr.: 0.1-1.2", "p" );
354 leg.AddEntry( gr3, "#mu-constr.: 0.1-1.4", "p" );
355 leg.AddEntry( gr4, "#mu-constr.: 0.2-1.1", "p" );
356 leg.AddEntry( gr5, "#mu-constr.: 0.2-1.3", "p" );
357 leg.AddEntry( gr6, "#mu-constr.: 0.2-1.5", "p" );
358 leg.AddEntry( gr7, "#mu-constr.: 0.3-1.2", "p" );
359 leg.AddEntry( gr8, "#mu-constr.: 0.3-1.4", "p" );
360 leg.AddEntry( gr9, "#mu-constr.: 0.3-1.6", "p" );
361
362 if not show_free_fit_res_34j:
363 leg.AddEntry( gr_1, "Constrained: 0.1-1.0", "p" );
364 leg.AddEntry( gr_2, "Constrained: 0.1-1.2", "p" );
365 leg.AddEntry( gr_3, "Constrained: 0.1-1.4", "p" );
366 leg.AddEntry( gr_4, "Constrained: 0.2-1.1", "p" );
367 leg.AddEntry( gr_5, "Constrained: 0.2-1.3", "p" );
368 leg.AddEntry( gr_6, "Constrained: 0.2-1.5", "p" );
369 leg.AddEntry( gr_7, "Constrained: 0.3-1.2", "p" );
370 leg.AddEntry( gr_8, "Constrained: 0.3-1.4", "p" );
371 leg.AddEntry( gr_9, "Constrained: 0.3-1.6", "p" );
372
373 leg.Draw();
374 if func == "gaus":
375 c2.SaveAs( "qcd_estimate_gaus.gif" );
376 elif( func == "pol3" ):
377 c2.SaveAs( "qcd_estimate_pol3.gif" );
378 elif( func == "landau" ):
379 c2.SaveAs( "qcd_estimate_landau.gif" );
380
381
382 h.GetYaxis().SetRangeUser( -1, 1 );
383 c2.SaveAs( Form( "qcd_estimate_%s_zoom_11.gif", func ) );
384
385 h.GetYaxis().SetRangeUser( -2, 2 );
386 c2.SaveAs( Form( "qcd_estimate_%s_zoom_22.gif", func ) );
387
388 h.GetYaxis().SetRangeUser( -3, 3 );
389 c2.SaveAs( Form( "qcd_estimate_%s_zoom_33.gif", func ) );
390
391 h.GetYaxis().SetRangeUser( -6, 6 );
392 c2.SaveAs( Form( "qcd_estimate_%s_zoom_66.gif", func ) );
393
394 h.GetYaxis().SetRangeUser( -8, 8 );
395 c2.SaveAs( Form( "qcd_estimate_%s_zoom_88.gif", func ) );
396
397
398
399 #------------------------------
400 # Plot 2: Average of 9 ranges
401 #------------------------------
402
403 if plot_average_est:
404 average_1j;
405 average_2j;
406 average_3j;
407 average_4mj;
408 average_3j_fix;
409 average_4mj_fix;
410
411 for i in range( 0, nrange ):
412 average_1j += y[i][0];
413 average_2j += y[i][1];
414 average_3j += y[i][2];
415 average_4mj += y[i][3];
416 average_3j_fix += yy34j[i][0];
417 average_4mj_fix += yy34j[i][1];
418
419 average_1j = average_1j / nrange;
420 average_2j = average_2j / nrange;
421 average_3j = average_3j / nrange;
422 average_4mj = average_4mj / nrange;
423 average_3j_fix = average_3j_fix / nrange;
424 average_4mj_fix = average_4mj_fix / nrange;
425 print "------------------------------------------" << endl;
426 print "average of 9 ranges, 1j: ", average_1j << endl;
427 print "average of 9 ranges, 2j: ", average_2j
428 print "average of 9 ranges, 3j: ", average_3j_fix
429 print "average of 9 ranges, 4mj: ", average_4mj_fix
430 print "average of 9 ranges, 3j (free): ", average_3j
431 print "average of 9 ranges, 4mj (free): ", average_4mj
432 print "------------------------------------------" << endl;
433
434 x12 = [1, 2]
435 av12 = [ average_1j, average_2j ]
436 av34 = [ average_3j, average_4mj ]
437 av34fix = [ average_3j_fix, average_4mj_fix]
438 gra12 = TGraph( 2, x12, av12 );
439 gra34 = TGraph( 2, xx, av34 );
440 gra34fix = TGraph( 2, xx, av34fix );
441
442
443 c3 = TCanvas( "c3", "Average QCD estimate", 610, 10, 600, 600 );
444 c3.cd();
445 c3.SetGrid( 1, 1 );
446 c3.SetLeftMargin( 0.12 );
447 c3.SetRightMargin( 0.35 );
448
449 hh = h.Clone();
450 hh.GetYaxis().SetRangeUser( -1, 1 );
451 hh.SetTitle( Form( "Average QCD estimate (%s)", func ) );
452 hh.Draw();
453
454 gra12.SetMarkerColor( kRed );
455 gra34.SetMarkerColor( kRed );
456 gra34fix.SetMarkerColor( kGreen + 3 );
457
458 gra12.SetMarkerStyle( 29 );
459 gra34.SetMarkerStyle( 30 );
460 gra34fix.SetMarkerStyle( 29 );
461
462 gra12.SetMarkerSize( 2.5 );
463 gra34.SetMarkerSize( 2.5 );
464 gra34fix.SetMarkerSize( 2.5 );
465
466 gra12.Draw( "P" );
467 gra34.Draw( "P" );
468 gra34fix.Draw( "P" );
469
470 leg2 = TLegend( 0.65, 0.7, 0.98, 0.9 );
471 leg2.SetFillColor( 0 );
472 leg2.AddEntry( gra12, "Free fits (1,2j)", "p" );
473 leg2.AddEntry( gra34, "Free fits (3,#geq4j)", "p" );
474 leg2.AddEntry( gra34fix, "Constrained fits (3,#geq4j)", "p" );
475 leg2.Draw();
476
477 c3.SaveAs( Form( "qcd_estimate_%s_average.gif", func ) );
478
479
480
481