Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
25void 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);
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/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,nullptr,"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}
#define f(i)
Definition RSha256.hxx:104
#define h(i)
Definition RSha256.hxx:106
int Int_t
Definition RtypesCore.h:45
float Float_t
Definition RtypesCore.h:57
double Double_t
Definition RtypesCore.h:59
@ kRed
Definition Rtypes.h:66
@ kBlue
Definition Rtypes.h:66
R__EXTERN TBenchmark * gBenchmark
Definition TBenchmark.h:59
#define gDirectory
Definition TDirectory.h:384
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
R__EXTERN TStyle * gStyle
Definition TStyle.h:433
virtual void Start(const char *name)
Starts Benchmark with the specified name.
virtual void Show(const char *name)
Stops Benchmark name and Prints results.
The Canvas class.
Definition TCanvas.h:23
virtual void ToggleEventStatus()
Toggle event statusbar.
Definition TCanvas.cxx:2425
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
Definition TFile.h:53
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition TFile.cxx:4067
static Bool_t SetCacheFileDir(std::string_view cacheDir, Bool_t operateDisconnected=kTRUE, Bool_t forceCacheread=kFALSE)
Sets the directory where to locally stage/cache remote files.
Definition TFile.cxx:4603
1-D histogram with a float per channel (see TH1 documentation)
Definition TH1.h:621
2D Histogram with Polygonal Bins
Definition TH2Poly.h:66
Book space in a file, create I/O buffers, to fill them, (un)compress them.
Definition TKey.h:28
virtual TObject * ReadObj()
To read a TObject* from the file.
Definition TKey.cxx:758
This class displays a legend box (TPaveText) containing several legend entries.
Definition TLegend.h:23
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition TMultiGraph.h:34
Mother of all ROOT objects.
Definition TObject.h:41
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition TObject.cxx:525
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:274
void SetGridx(Int_t value=1) override
Definition TPad.h:333
void SetGridy(Int_t value=1) override
Definition TPad.h:334
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
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:1636
void SetStatX(Float_t x=0)
Definition TStyle.h:394
void SetStatW(Float_t w=0.19)
Definition TStyle.h:396
void SetStatY(Float_t y=0)
Definition TStyle.h:395
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
leg
Definition legend1.C:34
#define H(x, y, z)
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Definition TMath.h:756
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
constexpr Double_t Pi()
Definition TMath.h:37
Double_t Tan(Double_t)
Returns the tangent of an angle of x radians.
Definition TMath.h:600
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123