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

# User Rev Content
1 dhidas 1.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