Logo ROOT   6.12/07
Reference Guide
tprofile2polyRealistic.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_hist
3 /// \notebook
4 /// Different charges depending on region
5 ///
6 /// \macro_image
7 /// \macro_code
8 ///
9 /// \author Filip Ilic
10 
11 #include <iostream>
12 #include <fstream>
13 #include <TProfile2Poly.h>
14 #include <TProfile2D.h>
15 #include <TCanvas.h>
16 #include <TRandom.h>
17 
18 using namespace std;
19 
20 void tprofile2polyRealistic(Int_t numEvents=100000)
21 {
22  int NUM_LS = 8;
23  TCanvas *c1 = new TCanvas("c1", "moving charge", 900, 400);
24  TCanvas *c2 = new TCanvas("c2", "Merge Individual moving charge plots", 800, 400);
25 
26  c1->Divide(NUM_LS, 3);
27  c2->Divide(3,1);
28 
29  // -------------------- Construct Reference plot bins ------------------------
30  auto new_avg = new TProfile2Poly();
31 
32  auto tot_avg_ls = new TProfile2Poly[NUM_LS];
33  auto det_avg_ls = new TProfile2Poly[NUM_LS];
34  auto det_err_ls = new TProfile2Poly[NUM_LS];
35  auto tot_merge = new TProfile2Poly();
36  auto det_avg_merge = new TProfile2Poly();
37  auto det_err_merge = new TProfile2Poly();
38 
39  float minx = -15;
40  float maxx = 15;
41  float miny = -15;
42  float maxy = 15;
43  float binsz = 0.5;
44 
45  for (float i = minx; i < maxx; i += binsz) {
46  for (float j = miny; j < maxy; j += binsz) {
47  tot_merge->AddBin(i, j, i + binsz, j + binsz);
48  for (int l=0; l<NUM_LS; ++l) {
49  tot_avg_ls[l].AddBin(i, j, i + binsz, j + binsz);
50  }
51  }
52  }
53 
54  // -------------------- Construct detector bins ------------------------
55  auto h2p = new TH2Poly();
56  auto tp2p = new TProfile2Poly();
57  ifstream infile;
58  TString dir = gROOT->GetTutorialDir();
59  dir.Append("/hist/data/tprofile2poly_tutorial.data");
60  infile.open(dir.Data());
61 
62  if (!infile) // Verify that the file was open successfully
63  {
64  std::cerr << dir.Data() << std::endl; // Report error
65  std::cerr << "Error code: " << strerror(errno) << std::endl; // Get some info as to why
66  return;
67  }
68  std::cout << " WE ARE AFTER LOADING DATA " << std::endl;
69 
70  vector<pair<Double_t, Double_t>> allCoords;
71  Double_t a, b;
72  while (infile >> a >> b) {
73  pair<Double_t, Double_t> coord(a, b);
74  allCoords.push_back(coord);
75  }
76 
77  if (allCoords.size() % 3 != 0) {
78  cout << "[ERROR] Bad file" << endl;
79  return;
80  }
81 
82  Double_t x[3], y[3];
83  for (Int_t i = 0; i < allCoords.size(); i += 3) {
84  x[0] = allCoords[i + 0].first;
85  y[0] = allCoords[i + 0].second;
86  x[1] = allCoords[i + 1].first;
87  y[1] = allCoords[i + 1].second;
88  x[2] = allCoords[i + 2].first;
89  y[2] = allCoords[i + 2].second;
90 
91  det_avg_merge->AddBin(3, x, y);
92  det_err_merge->AddBin(3, x, y);
93 
94  for (int l=0; l<NUM_LS; ++l) {
95  det_avg_ls[l].AddBin(3, x, y);
96  }
97  }
98 
99  std::cout << " WE ARE AFTER ADDING BINS " << std::endl;
100 
101  // -------------------- Simulate particles ------------------------
102  TRandom ran;
103 
104  // moving error
105  Double_t xoffset1 = 0;
106  Double_t yoffset1 = 0;
107  Double_t xoffset2 = 0;
108  Double_t yoffset2 = 0;
109 
110  for (int i = 0; i <= NUM_LS-1; ++i) { // LumiSection
111  std::cout << "[In Progress] LumiSection " << i << std::endl;
112  for (int j = 0; j < numEvents; ++j) { // Events
113  Double_t r1 = ran.Gaus(0, 10);
114  Double_t r2 = ran.Gaus(0, 8);
115  Double_t rok = ran.Gaus(10, 1);
116  Double_t rbad1 = ran.Gaus(8, 5);
117  Double_t rbad2 = ran.Gaus(-8, 5);
118 
119  Double_t val = rok;
120 
121  xoffset1 += 0.00002;
122  yoffset1 += 0.00002;
123 
124  xoffset2 += 0.00003;
125  yoffset2 += 0.00004;
126 
127  if (r2 > 3. - yoffset1 && r2 < 8. - yoffset1 &&
128  r1 > 1. + xoffset1 && r1 < 5. + xoffset1 ) {
129  val -= rbad1;
130  }
131 
132  if (r2 > -10 + yoffset2 && r2 < -8 + yoffset2 &&
133  r1 > -6 + xoffset2 && r1 < 8 + xoffset2 ) {
134  val -= rbad2;
135  }
136 
137  tot_avg_ls[i].Fill(r1, r2, val);
138  det_avg_ls[i].Fill(r1, r2, val);
139  }
140 
141  std::string title;
142 
143 
144  c1->cd(i+1);
145  title = "Global View: Avg in LS " + to_string(i);
146  tot_avg_ls[i].SetTitle(title.c_str());
147  tot_avg_ls[i].SetStats(false);
148  tot_avg_ls[i].Draw("COLZ");
149  c1->Update();
150 
151  c1->cd((i+1)+NUM_LS);
152  title = "Detector View: Avg in LS " + to_string(i);
153  det_avg_ls[i].SetTitle(title.c_str());
154  det_avg_ls[i].SetStats(false);
155  det_avg_ls[i].Draw("COLZ");
156  c1->Update();
157 
158 
159  c1->cd((i+1)+(NUM_LS*2));
160  title = "Detector View: Error in LS " + to_string(i);
161  det_avg_ls[i].SetTitle(title.c_str());
162  det_avg_ls[i].SetStats(false);
163  det_avg_ls[i].SetContentToError();
164  det_avg_ls[i].Draw("COLZ");
165  c1->Update();
166  }
167 
168  std::vector<TProfile2Poly*> tot_avg_v;
169  std::vector<TProfile2Poly*> det_avg_v;
170  for (Int_t t=0; t<NUM_LS; t++){
171  tot_avg_v.push_back(&tot_avg_ls[t]);
172  det_avg_v.push_back(&det_avg_ls[t]);
173  }
174 
175  std::cout << "[In Progress] Merging" << std::endl;
176 
177  std::string title;
178 
179  tot_merge->Merge(tot_avg_v);
180  c2->cd(1);
181  title = "Total average merge";
182  tot_merge->SetTitle(title.c_str());
183  tot_merge->Draw("COLZ");
184 
185  det_avg_merge->Merge(det_avg_v);
186  c2->cd(2);
187  title = "Detector average merge";
188  det_avg_merge->SetTitle(title.c_str());
189  det_avg_merge->SetContentToAverage(); // implicit
190  det_avg_merge->Draw("COLZ");
191 
192  det_err_merge->Merge(det_avg_v);
193  c2->cd(3);
194  title = "Detector error merge";
195  det_err_merge->SetTitle(title.c_str());
196  det_err_merge->SetContentToError();
197  det_err_merge->Draw("COLZ");
198 }
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:256
return c1
Definition: legend1.C:41
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
#define gROOT
Definition: TROOT.h:402
int Int_t
Definition: RtypesCore.h:41
STL namespace.
Double_t x[n]
Definition: legend1.C:17
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
auto * a
Definition: textangle.C:12
2D Profile Histogram with Polygonal Bins.
Definition: TProfile2Poly.h:57
virtual void SetTitle(const char *title="")=0
The Canvas class.
Definition: TCanvas.h:31
return c2
Definition: legend2.C:14
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
auto * l
Definition: textangle.C:4
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2248
2D Histogram with Polygonal Bins
Definition: TH2Poly.h:66