Logo ROOT   6.18/05
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
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.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
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}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define f(i)
Definition: RSha256.hxx:104
#define h(i)
Definition: RSha256.hxx:106
#define R(a, b, c, d, e, f, g, h, i)
Definition: RSha256.hxx:110
int Int_t
Definition: RtypesCore.h:41
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
@ kRed
Definition: Rtypes.h:64
@ kBlue
Definition: Rtypes.h:64
R__EXTERN TBenchmark * gBenchmark
Definition: TBenchmark.h:59
#define gDirectory
Definition: TDirectory.h:218
R__EXTERN TStyle * gStyle
Definition: TStyle.h:406
virtual void SetNdivisions(Int_t n=510, Bool_t optim=kTRUE)
Set the number of divisions for this axis.
Definition: TAttAxis.cxx:229
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
virtual void Start(const char *name)
Starts Benchmark with the specified name.
Definition: TBenchmark.cxx:174
virtual void Show(const char *name)
Stops Benchmark name and Prints results.
Definition: TBenchmark.cxx:157
The Canvas class.
Definition: TCanvas.h:31
virtual void ToggleEventStatus()
Toggle event statusbar.
Definition: TCanvas.cxx:2228
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:48
static Bool_t SetCacheFileDir(ROOT::Internal::TStringView cacheDir, Bool_t operateDisconnected=kTRUE, Bool_t forceCacheread=kFALSE)
Definition: TFile.h:316
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseGeneralPurpose, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:3980
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:571
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9519
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8476
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
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:7964
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:8619
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3258
TAxis * GetYaxis()
Definition: TH1.h:317
virtual void SetContour(Int_t nlevels, const Double_t *levels=0)
Set the number and values of contour levels.
Definition: TH1.cxx:7905
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7413
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6218
2D Histogram with Polygonal Bins
Definition: TH2Poly.h:66
Int_t GetNumberOfBins() const
Definition: TH2Poly.h:101
virtual Int_t Fill(Double_t x, Double_t y)
Increment the bin containing (x,y) by 1.
Definition: TH2Poly.cxx:616
void ChangePartition(Int_t n, Int_t m)
Changes the number of partition cells in the histogram.
Definition: TH2Poly.cxx:467
const char * GetBinName(Int_t bin) const
Returns the bin name.
Definition: TH2Poly.cxx:845
virtual Double_t GetBinContent(Int_t bin) const
Returns the content of the input bin For the overflow/underflow/sea bins:
Definition: TH2Poly.cxx:788
virtual Int_t AddBin(TObject *poly)
Adds a new bin to the histogram.
Definition: TH2Poly.cxx:222
virtual Double_t GetBinError(Int_t bin) const
Returns the value of error associated to bin number bin.
Definition: TH2Poly.cxx:803
virtual 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:1303
Double_t GetMaximum() const
Returns the maximum value of the histogram.
Definition: TH2Poly.cxx:865
Book space in a file, create I/O buffers, to fill them, (un)compress them.
Definition: TKey.h:24
virtual TObject * ReadObj()
To read a TObject* from the file.
Definition: TKey.cxx:722
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:35
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
Mother of all ROOT objects.
Definition: TObject.h:37
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:443
virtual void SetGridx(Int_t value=1)
Definition: TPad.h:329
virtual void SetGridy(Int_t value=1)
Definition: TPad.h:330
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:1444
void SetStatX(Float_t x=0)
Definition: TStyle.h:374
void SetStatW(Float_t w=0.19)
Definition: TStyle.h:376
void SetStatY(Float_t y=0)
Definition: TStyle.h:375
return c1
Definition: legend1.C:41
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
leg
Definition: legend1.C:34
#define H(x, y, z)
static constexpr double mg
Double_t Log(Double_t x)
Definition: TMath.h:748
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
constexpr Double_t Pi()
Definition: TMath.h:38
Double_t Tan(Double_t)
Definition: TMath.h:633
Short_t Abs(Short_t d)
Definition: TMathBase.h:120