Logo ROOT   6.12/07
Reference Guide
th2polyEurope.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_hist
3 /// \notebook -js
4 /// This tutorial illustrates how to create an histogram with polygonal
5 /// bins (TH2Poly), fill it and draw it. The initial data are stored
6 /// in TMultiGraphs. They represent the european countries.
7 /// The histogram filling is done according to a Mercator projection,
8 /// therefore the bin contains should be proportional to the real surface
9 /// of the countries.
10 ///
11 /// The initial data have been downloaded from: http://www.maproom.psu.edu/dcw/
12 /// This database was developed in 1991/1992 and national boundaries reflect
13 /// political reality as of that time.
14 ///
15 /// The script is shooting npoints (script argument) randomly over the Europe area.
16 /// The number of points inside the countries should be proportional to the country surface
17 /// The estimated surface is compared to the surfaces taken from wikipedia.
18 ///
19 /// \macro_image
20 /// \macro_output
21 /// \macro_code
22 ///
23 /// \author Olivier Couet
24 
25 void th2polyEurope(Int_t npoints=500000)
26 {
27  Int_t i,j;
28  Double_t lon1 = -25;
29  Double_t lon2 = 35;
30  Double_t lat1 = 34;
31  Double_t lat2 = 72;
32  Double_t R = (lat2-lat1)/(lon2-lon1);
33  Int_t W = 800;
34  Int_t H = (Int_t)(R*800);
35  gStyle->SetStatX(0.28);
36  gStyle->SetStatY(0.45);
37  gStyle->SetStatW(0.15);
38 
39  // Canvas used to draw TH2Poly (the map)
40  TCanvas *ce = new TCanvas("ce", "ce",0,0,W,H);
41  ce->ToggleEventStatus();
42  ce->SetGridx();
43  ce->SetGridy();
44 
45  // Real surfaces taken from Wikipedia.
46  const Int_t nx = 36;
47  // see http://en.wikipedia.org/wiki/Area_and_population_of_European_countries
48  const char *countries[nx] = {
49  "france", "spain", "sweden", "germany", "finland",
50  "norway", "poland", "italy", "yugoslavia", "united_kingdom",
51  "romania", "belarus","greece", "czechoslovakia","bulgaria",
52  "iceland", "hungary","portugal","austria", "ireland",
53  "lithuania", "latvia", "estonia", "denmark", "netherlands",
54  "switzerland","moldova","belgium", "albania", "cyprus",
55  "luxembourg", "andorra","malta", "liechtenstein", "san_marino", "monaco" };
56  Float_t surfaces[nx] = {
57  547030, 505580, 449964, 357021, 338145,
58  324220, 312685, 301230, 255438, 244820,
59  237500, 207600, 131940, 127711, 110910,
60  103000, 93030, 89242, 83870, 70280,
61  65200, 64589, 45226, 43094, 41526,
62  41290, 33843, 30528, 28748, 9250,
63  2586, 468, 316, 160, 61, 2};
64 
65  TH1F *h = new TH1F("h","Countries surfaces (in km^{2})",3,0,3);
66  for (i=0; i<nx; i++) h->Fill(countries[i], surfaces[i]);
67  h->LabelsDeflate();
68 
70  TFile *f;
71  f = TFile::Open("http://root.cern.ch/files/europe.root","cacheread");
72 
73  if (!f) {
74  printf("Cannot access europe.root. Is internet working ?\n");
75  return;
76  }
77 
78  TH2Poly *p = new TH2Poly(
79  "Europe",
80  "Europe (bin contents are normalized to the surfaces in km^{2})",
81  lon1,lon2,lat1,lat2);
82  p->GetXaxis()->SetNdivisions(520);
83  p->GetXaxis()->SetTitle("longitude");
84  p->GetYaxis()->SetTitle("latitude");
85 
86  p->SetContour(100);
87 
88  TMultiGraph *mg;
89  TKey *key;
90  TIter nextkey(gDirectory->GetListOfKeys());
91  while ((key = (TKey*)nextkey())) {
92  TObject *obj = key->ReadObj();
93  if (obj->InheritsFrom("TMultiGraph")) {
94  mg = (TMultiGraph*)obj;
95  p->AddBin(mg);
96  }
97  }
98 
99  TRandom r;
100  Double_t longitude, latitude;
101  Double_t x, y, pi4 = TMath::Pi()/4, alpha = TMath::Pi()/360;
102 
103  gBenchmark->Start("Partitioning");
104  p->ChangePartition(100, 100);
105  gBenchmark->Show("Partitioning");
106 
107  // Fill TH2Poly according to a Mercator projection.
108  gBenchmark->Start("Filling");
109  for (i=0; i<npoints; i++) {
110  longitude = r.Uniform(lon1,lon2);
111  latitude = r.Uniform(lat1,lat2);
112  x = longitude;
113  y = 38*TMath::Log(TMath::Tan(pi4+alpha*latitude));
114  p->Fill(x,y);
115  }
116  gBenchmark->Show("Filling");
117 
118  Int_t nbins = p->GetNumberOfBins();
119  Double_t maximum = p->GetMaximum();
120 
121 
122  // h2 contains the surfaces computed from TH2Poly.
123  TH1F *h2 = (TH1F *)h->Clone("h2");
124  h2->Reset();
125  for (j=0; j<nx; j++) {
126  for (i=0; i<nbins; i++) {
127  if (strstr(countries[j],p->GetBinName(i+1))) {
128  h2->Fill(countries[j],p->GetBinContent(i+1));
129  h2->SetBinError(j, p->GetBinError(i+1));
130  }
131  }
132  }
133 
134  // Normalize the TH2Poly bin contents to the real surfaces.
135  Double_t scale = surfaces[0]/maximum;
136  for (i=0; i<nbins; i++) p->SetBinContent(i+1, scale*p->GetBinContent(i+1));
137 
138  gStyle->SetOptStat(1111);
139  p->Draw("COLZ");
140 
141  TCanvas *c1 = new TCanvas("c1", "c1",W+10,0,W-20,H);
142  c1->SetRightMargin(0.047);
143 
144  scale = h->GetMaximum()/h2->GetMaximum();
145 
146  h->SetStats(0);
147  h->SetLineColor(kRed-3);
148  h->SetLineWidth(2);
149  h->SetMarkerStyle(20);
150  h->SetMarkerColor(kBlue);
151  h->SetMarkerSize(0.8);
152  h->Draw("LP");
153  h->GetXaxis()->SetLabelFont(42);
154  h->GetXaxis()->SetLabelSize(0.03);
155  h->GetYaxis()->SetLabelFont(42);
156 
157  h2->Scale(scale);
158  Double_t scale2=TMath::Sqrt(scale);
159  for (i=0; i<nx; i++) h2->SetBinError(i+1, scale2*h2->GetBinError(i+1));
160  h2->Draw("E SAME");
161  h2->SetMarkerStyle(20);
162  h2->SetMarkerSize(0.8);
163 
164  TLegend *leg = new TLegend(0.5,0.67,0.92,0.8,NULL,"NDC");
165  leg->SetTextFont(42);
166  leg->SetTextSize(0.025);
167  leg->AddEntry(h,"Real countries surfaces from Wikipedia (in km^{2})","lp");
168  leg->AddEntry(h2,"Countries surfaces from TH2Poly (with errors)","lp");
169  leg->Draw();
170  leg->Draw();
171 
172  Double_t wikiSum = h->Integral();
173  Double_t polySum = h2->Integral();
174  Double_t error = TMath::Abs(wikiSum-polySum)/wikiSum;
175  printf("THPoly Europe surface estimation error wrt wikipedia = %f per cent when using %d points\n",100*error,npoints);
176 }
virtual void SetGridx(Int_t value=1)
Definition: TPad.h:328
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6063
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3244
virtual Double_t GetMaximum(Double_t maxval=FLT_MAX) const
Return maximum value smaller than maxval of bins in the range, unless the value has been overridden b...
Definition: TH1.cxx:7807
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
virtual void ToggleEventStatus()
Toggle event statusbar.
Definition: TCanvas.cxx:2190
Double_t Log(Double_t x)
Definition: TMath.h:648
virtual Double_t GetBinError(Int_t bin) const
Returns the value of error associated to bin number bin.
Definition: TH2Poly.cxx:775
float Float_t
Definition: RtypesCore.h:53
return c1
Definition: legend1.C:41
Definition: Rtypes.h:59
void ChangePartition(Int_t n, Int_t m)
Changes the number of partition cells in the histogram.
Definition: TH2Poly.cxx:457
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition: TH1.cxx:7748
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:452
THist< 1, float, THistStatContent, THistStatUncertainty > TH1F
Definition: THist.hxx:285
TH1 * h
Definition: legend2.C:5
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:46
virtual void SetNdivisions(Int_t n=510, Bool_t optim=kTRUE)
Set the number of divisions for this axis.
Definition: TAttAxis.cxx:229
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:35
static Bool_t SetCacheFileDir(ROOT::Internal::TStringView cacheDir, Bool_t operateDisconnected=kTRUE, Bool_t forceCacheread=kFALSE)
Definition: TFile.h:303
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7256
#define H(x, y, z)
virtual void Show(const char *name)
Stops Benchmark name and Prints results.
Definition: TBenchmark.cxx:157
void SetStatX(Float_t x=0)
Definition: TStyle.h:370
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:567
int Int_t
Definition: RtypesCore.h:41
void SetStatY(Float_t y=0)
Definition: TStyle.h:371
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
static constexpr double mg
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9361
virtual void SetLabelFont(Style_t font=62)
Set labels&#39; font.
Definition: TAttAxis.cxx:183
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=1, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3950
virtual void SetTextFont(Font_t tfont=62)
Set the text font.
Definition: TAttText.h:45
virtual void Start(const char *name)
Starts Benchmark with the specified name.
Definition: TBenchmark.cxx:174
Double_t x[n]
Definition: legend1.C:17
const char * GetBinName(Int_t bin) const
Returns the bin name.
Definition: TH2Poly.cxx:791
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual Int_t AddBin(TObject *poly)
Adds a new bin to the histogram.
Definition: TH2Poly.cxx:215
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
Book space in a file, create I/O buffers, to fill them, (un)compress them.
Definition: TKey.h:24
constexpr Double_t Pi()
Definition: TMath.h:40
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:8461
virtual Double_t GetBinContent(Int_t bin) const
Returns the content of the input bin For the overflow/underflow/sea bins: -1 | -2 | -3 ---+----+---- ...
Definition: TH2Poly.cxx:762
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
void SetBinContent(Int_t bin, Double_t content)
Sets the contents of the input bin to the input content Negative values between -1 and -9 are for the...
Definition: TH2Poly.cxx:1248
R__EXTERN TBenchmark * gBenchmark
Definition: TBenchmark.h:59
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
virtual void LabelsDeflate(Option_t *axis="X")
Reduce the number of bins for the axis passed in the option to the number of bins having a label...
Definition: TH1.cxx:4928
Int_t GetNumberOfBins() const
Definition: TH2Poly.h:110
void SetStatW(Float_t w=0.19)
Definition: TStyle.h:372
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
TAxis * GetYaxis()
Definition: TH1.h:316
virtual void SetLabelSize(Float_t size=0.04)
Set size of axis labels The size is expressed in per cent of the pad width.
Definition: TAttAxis.cxx:204
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
The Canvas class.
Definition: TCanvas.h:31
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:359
virtual void SetGridy(Int_t value=1)
Definition: TPad.h:329
leg
Definition: legend1.C:34
Double_t y[n]
Definition: legend1.C:17
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:627
virtual void SetRightMargin(Float_t rightmargin)
Set Pad right margin in fraction of the pad width.
Definition: TAttPad.cxx:120
Mother of all ROOT objects.
Definition: TObject.h:37
virtual TObject * ReadObj()
To read a TObject* from the file.
Definition: TKey.cxx:722
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2662
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1266
#define gDirectory
Definition: TDirectory.h:213
Definition: Rtypes.h:59
Double_t Sqrt(Double_t x)
Definition: TMath.h:590
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition: TAttText.h:46
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
constexpr Double_t R()
Definition: TMath.h:213
Double_t GetMaximum() const
Returns the maximum value of the histogram.
Definition: TH2Poly.cxx:811
virtual void SetStats(Bool_t stats=kTRUE)
Set statistics option on/off.
Definition: TH1.cxx:8247
Double_t Tan(Double_t)
Definition: TMath.h:553
virtual Int_t Fill(Double_t x, Double_t y)
Increment the bin containing (x,y) by 1.
Definition: TH2Poly.cxx:605
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:315
2D Histogram with Polygonal Bins
Definition: TH2Poly.h:66
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8319