Belle II Software  release-05-01-25
ARICHChannelHist.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2016 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Luka Santelj *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <arich/utility/ARICHChannelHist.h>
12 #include <framework/logging/Logger.h>
13 #include <math.h>
14 #include <iostream>
15 #include <algorithm>
16 #include <TGraph.h>
17 #include <TVector2.h>
18 #include <TList.h>
19 
20 using namespace std;
21 using namespace Belle2;
22 
23 ARICHChannelHist::ARICHChannelHist(const char* name, const char* title, int type,
24  const std::vector<unsigned>& moduleIDs) : TH2Poly()
25 {
26 
27  m_type = type;
28  SetName(name);
29  SetTitle(title);
30  m_hapd2binMap.assign(420, 0);
31 
32  // positions of HAPDs and channel mapping (avoid using DB classes...)
33  double rs[7] = {57.35, 65.81, 74.37, 82.868, 91.305, 99.794, 108.185};
34  unsigned nhapds[7] = {42, 48, 54, 60, 66, 72, 78};
35  unsigned chmap[144] = {88, 86, 96, 87, 75, 72, 97, 108, 73, 74, 98, 109, 84, 85, 120, 110, 76, 77, 132, 121, 136, 124, 99, 133, 125, 113, 122, 111, 101, 137, 134, 123, 89, 112, 100, 135, 52, 28, 3, 40, 41, 5, 15, 2, 17, 29, 27, 14, 4, 16, 1, 26, 53, 65, 0, 13, 48, 49, 39, 12, 61, 62, 25, 38, 63, 60, 24, 37, 64, 50, 51, 36, 91, 115, 140, 103, 102, 138, 128, 141, 126, 114, 116, 129, 139, 127, 142, 117, 90, 78, 143, 130, 95, 94, 104, 131, 82, 81, 118, 105, 80, 83, 119, 106, 79, 93, 92, 107, 55, 57, 47, 56, 68, 71, 46, 35, 70, 69, 45, 34, 59, 58, 23, 33, 67, 66, 11, 22, 7, 19, 44, 10, 18, 30, 21, 32, 42, 6, 9, 20, 54, 31, 43, 8};
36  double chns[12] = { -2.88, -2.37, -1.86, -1.35, -0.84, -0.33, 0.33, 0.84, 1.35, 1.86, 2.37, 2.88};
37 
38  float size = 0.5 / 2. - 0.01;
39  if (m_type == 1) size = 7.0 / 2. - 0.5;
40  if (m_type == 2) size = 3.3 / 2.;
41 
42  float X[5], Y[5], globX[5], globY[5];
43  X[0] = -size; Y[0] = -size;
44  X[1] = size; Y[1] = -size;
45  X[2] = size; Y[2] = size;
46  X[3] = -size; Y[3] = size;
47  X[4] = -size; Y[4] = -size;
48 
49  int nhapd = 1;
50  unsigned iring = 0;
51  unsigned ihapd = 0;
52 
53  std::vector<unsigned> ids;
54  if (moduleIDs.size() > 0) ids = moduleIDs;
55  else {
56  for (int hapdID = 1; hapdID < 421; hapdID++) {
57  ids.push_back(hapdID);
58  }
59  }
60 
61  // HAPD bins
62  if (m_type == 1) {
63  for (int hapdID = 1; hapdID < 421; hapdID++) {
64  //for (unsigned hapdID : ids) {
65  //m_hapd2binMap[hapdID - 1] = nhapd;
66  //nhapd++;
67  float r = rs[iring];
68  float dphi = 2.*M_PI / nhapds[iring];
69  float fi = dphi / 2. + ihapd * dphi;
70  TVector2 centerPos(r * cos(fi), r * sin(fi));
71  for (int i = 0; i < 5; i++) {
72  float rotX = X[i] * cos(fi) - Y[i] * sin(fi);
73  float rotY = X[i] * sin(fi) + Y[i] * cos(fi);
74  globX[i] = rotX + centerPos.X();
75  globY[i] = rotY + centerPos.Y();
76  }
77  if (std::find(ids.begin(), ids.end(), hapdID) != ids.end()) {
78  m_hapd2binMap[hapdID - 1] = nhapd;
79  nhapd++;
80  TGraph* mybox = new TGraph(5, globX, globY);
81  mybox->SetName((to_string(hapdID)).c_str());
82  AddBin(mybox);
83  }
84  ihapd++;
85  if (ihapd == nhapds[iring]) { iring++; ihapd = 0;}
86  }
87 
88  } else if (m_type == 0) {
89  for (int hapdID = 1; hapdID < 421; hapdID++) {
90  //for (unsigned hapdID : ids) {
91  // m_hapd2binMap[hapdID - 1] = nhapd;
92  // nhapd++;
93  float dphi = 2.*M_PI / nhapds[iring];
94  float fi = dphi / 2. + ihapd * dphi;
95  float r = rs[iring];
96  for (int chID = 0; chID < 144; chID++) {
97 
98  unsigned chX = chmap[chID] % 12;
99  unsigned chY = chmap[chID] / 12;
100  TVector2 hapdPos(r * cos(fi), r * sin(fi));
101  TVector2 locPos(chns[chX], chns[chY]);
102  TVector2 centerPos = hapdPos + locPos.Rotate(fi);
103 
104  for (int i = 0; i < 5; i++) {
105  float rotX = X[i] * cos(fi) - Y[i] * sin(fi);
106  float rotY = X[i] * sin(fi) + Y[i] * cos(fi);
107  globX[i] = rotX + centerPos.X();
108  globY[i] = rotY + centerPos.Y();
109  }
110  if (std::find(ids.begin(), ids.end(), hapdID) != ids.end()) {
111  m_hapd2binMap[hapdID - 1] = nhapd;
112  if (chID == 143) nhapd++;
113  TGraph* mybox = new TGraph(5, globX, globY);
114  mybox->SetName((to_string(hapdID)).c_str());
115  AddBin(mybox);
116  }
117  }
118  ihapd++;
119  if (ihapd == nhapds[iring]) { iring++; ihapd = 0;}
120  }
121  } else if (m_type == 2) {
122  size += 0.2;
123  for (int hapdID = 1; hapdID < 421; hapdID++) {
124  float dphi = 2.*M_PI / nhapds[iring];
125  float fi = dphi / 2. + ihapd * dphi;
126  float r = rs[iring];
127  TVector2 hapdPos(r * cos(fi), r * sin(fi));
128  for (int chipID = 0; chipID < 4; chipID++) {
129  TVector2 locPos(-size + (chipID / 2)*size * 2, size - (chipID % 2)*size * 2);
130  TVector2 centerPos = hapdPos + locPos.Rotate(fi);
131 
132  for (int i = 0; i < 5; i++) {
133  float rotX = X[i] * cos(fi) - Y[i] * sin(fi);
134  float rotY = X[i] * sin(fi) + Y[i] * cos(fi);
135  globX[i] = rotX + centerPos.X();
136  globY[i] = rotY + centerPos.Y();
137  }
138 
139  if (std::find(ids.begin(), ids.end(), hapdID) != ids.end()) {
140  m_hapd2binMap[hapdID - 1] = nhapd;
141  if (chipID == 3) nhapd++;
142  TGraph* mybox = new TGraph(5, globX, globY);
143  mybox->SetName((to_string(hapdID)).c_str());
144  AddBin(mybox);
145  }
146  }
147  ihapd++;
148  if (ihapd == nhapds[iring]) { iring++; ihapd = 0;}
149  }
150  } else std::cout << "Invalid histogram type! use 0 for channel bins or 1 for HAPD bins" << std::endl;
151  SetOption("colz");
152  SetStats(0);
153  GetXaxis()->SetLimits(-115., 115.);
154  GetYaxis()->SetLimits(-115., 115.);
155 }
156 
157 void ARICHChannelHist::Draw(Option_t* option)
158 {
159  TH2Poly::Draw(option);
160  double rlin = 40;
161  double rlout = 113;
162  for (int isec = 0; isec < 6; isec++) {
163  double x1 = rlin * cos(M_PI / 3.*isec);
164  double x2 = rlout * cos(M_PI / 3.*isec);
165  double y1 = rlin * sin(M_PI / 3.*isec);
166  double y2 = rlout * sin(M_PI / 3.*isec);
167  lines[isec] = TLine(x1, y1, x2, y2);
168  lines[isec].Draw();
169  x1 = rlin * cos(M_PI / 3.*isec + M_PI / 6.);
170  y1 = rlin * sin(M_PI / 3.*isec + M_PI / 6.);
171  labels[isec] = TText(x1, y1, TString::Format("S-%d", isec + 1));
172  labels[isec].SetTextAlign(22);
173  labels[isec].SetTextSize(0.03);
174  labels[isec].Draw();
175  }
176 }
177 
178 
179 void ARICHChannelHist::fillBin(unsigned hapdID, unsigned chID, double weight)
180 {
181  unsigned chIndex = 0;
182  if (m_type == 0) chIndex = (m_hapd2binMap[hapdID - 1] - 1) * 144 + chID + 1;
183  if (m_type == 2) chIndex = (m_hapd2binMap[hapdID - 1] - 1) * 4 + chID + 1;
184  SetBinContent(chIndex, GetBinContent(chIndex) + weight);
185 }
186 
187 void ARICHChannelHist::setBinContent(unsigned hapdID, unsigned chID, double value)
188 {
189  unsigned chIndex = 0;
190  if (m_type == 0) chIndex = (m_hapd2binMap[hapdID - 1] - 1) * 144 + chID + 1;
191  if (m_type == 2) chIndex = (m_hapd2binMap[hapdID - 1] - 1) * 4 + chID + 1;
192  SetBinContent(chIndex, value);
193 }
194 
195 void ARICHChannelHist::setBinContent(unsigned hapdID, double value)
196 {
197  SetBinContent(m_hapd2binMap[hapdID - 1], value);
198 }
199 
200 void ARICHChannelHist::fillBin(unsigned hapdID, double weight)
201 {
202  SetBinContent(m_hapd2binMap[hapdID - 1], GetBinContent(hapdID) + weight);
203 }
204 
205 void ARICHChannelHist::fillFromTH1(TH1* hist)
206 {
207 
208  int nbins = hist->GetNbinsX();
209  if (m_type == 1) {
210  if (nbins < 420) { B2ERROR("Number of bins in histogram small than number of ChannelHist bins!"); return;}
211  if (nbins == 420) for (int i = 0; i < nbins; i++) setBinContent(i + 1, hist->GetBinContent(i + 1));
212  if (nbins == 420 * 4) {
213  for (int i = 0; i < 420; i++) {
214  for (int j = 0; j < 4; j++) fillBin(i + 1, hist->GetBinContent(i * 4 + j + 1));
215  }
216  }
217  if (nbins == 420 * 144) {
218  for (int i = 0; i < 420; i++) {
219  for (int j = 0; j < 144; j++) fillBin(i + 1, hist->GetBinContent(i * 144 + j + 1));
220  }
221  }
222  } else if (m_type == 0) {
223  if (nbins < 420 * 144) { B2ERROR("Number of bins in histogram small than number of ChannelHist bins!"); return;}
224  for (int i = 0; i < 420; i++) {
225  for (int j = 0; j < 144; j++) setBinContent(i + 1, j, hist->GetBinContent(i * 144 + j + 1));
226  }
227  } else if (m_type == 2) {
228  if (nbins < 420 * 4) { B2ERROR("Number of bins in histogram small than number of ChannelHist bins!"); return;}
229  if (nbins == 420 * 4) for (int i = 0; i < nbins; i++) setBinContent(i / 4 + 1, i % 4, hist->GetBinContent(i + 1));
230  if (nbins == 420 * 144) {
231  for (int i = 0; i < 420; i++) {
232  for (int j = 0; j < 144; j++) fillBin(i + 1, j / 36, hist->GetBinContent(i * 144 + j + 1));
233  }
234  }
235  } else return;
236 }
237 
238 void ARICHChannelHist::setPoly(TH2Poly* poly)
239 {
240 
241  if (poly->GetNumberOfBins() == 0) {
242  for (const auto && bin : *fBins) {
243  poly->AddBin((TGraph*)((TH2PolyBin*)bin)->GetPolygon());
244  }
245  }
246  if (poly->GetNumberOfBins() != GetNumberOfBins()) {std::cout << "Mismatch between number of bins in TH2Poly and ARICHChannelHist" << std::endl; return;}
247 
248  double max = poly->GetMaximum();
249  for (int i = 1; i < GetNumberOfBins() + 1; i++) {
250  poly->SetBinContent(i, GetBinContent(i) > max ? max : GetBinContent(i));
251  }
252  return;
253 }
254 
Belle2::ARICHChannelHist::fillBin
void fillBin(unsigned hapdID, unsigned chID, double weight=1.)
Add entry to bin corresponding to hapd hapdID and channel chID.
Definition: ARICHChannelHist.cc:179
Belle2::ARICHChannelHist::Draw
void Draw(Option_t *option="") override
Draw the histogram.
Definition: ARICHChannelHist.cc:157
Belle2::ARICHChannelHist::m_type
int m_type
histogram type
Definition: ARICHChannelHist.h:109
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ARICHChannelHist::setPoly
void setPoly(TH2Poly *poly)
Fill pure TH2Poly from ARICHChannelHist, makes bins and fills content.
Definition: ARICHChannelHist.cc:238
Belle2::ARICHChannelHist::m_hapd2binMap
std::vector< unsigned > m_hapd2binMap
map of bins
Definition: ARICHChannelHist.h:110
Belle2::ARICHChannelHist::fillFromTH1
void fillFromTH1(TH1 *hist)
Fill the channelHist from the histogram Type 0 channelHist has to be filled with 420*144bin TH1 (each...
Definition: ARICHChannelHist.cc:205
Belle2::ARICHChannelHist::setBinContent
void setBinContent(unsigned hapdID, unsigned chID, double value)
Set content of bin corresponding to hapd hapdID and channel chID.
Definition: ARICHChannelHist.cc:187