ATTPCROOT  0.3.0-alpha
A ROOT-based framework for analyzing data from active target detectors
AtSpecMATMap.cxx
Go to the documentation of this file.
1 #include "AtSpecMATMap.h"
2 
3 #include <FairLogger.h>
4 
5 #include <Math/Point2D.h>
6 #include <Rtypes.h>
7 #include <TAxis.h>
8 #include <TH2Poly.h>
9 
10 #include <boost/multi_array/base.hpp>
11 #include <boost/multi_array/extent_gen.hpp>
12 #include <boost/multi_array/multi_array_ref.hpp>
13 #include <boost/multi_array/subarray.hpp>
14 
15 #include <algorithm>
16 #include <cmath>
17 #include <fstream> // IWYU pragma: keep
18 #include <iostream>
19 #include <vector>
20 
21 constexpr auto cRED = "\033[1;31m";
22 constexpr auto cYELLOW = "\033[1;33m";
23 constexpr auto cNORMAL = "\033[0m";
24 
25 using XYPoint = ROOT::Math::XYPoint;
26 
27 AtSpecMATMap::AtSpecMATMap(Int_t numPads) : AtMap()
28 {
29  AtPadCoord.resize(boost::extents[numPads][4][2]);
30  std::fill(AtPadCoord.data(), AtPadCoord.data() + AtPadCoord.num_elements(), 0);
31  std::cout << " SpecMAT Map initialized " << std::endl;
32  std::cout << " SpecMAT Pad Coordinates container initialized " << std::endl;
33  fNumberPads = numPads;
34 }
35 
36 AtSpecMATMap::~AtSpecMATMap() = default;
37 
38 void AtSpecMATMap::Dump()
39 {
40 
41  std::ofstream coordmap;
42  coordmap.open("coordmap_specmat.txt");
43  coordmap << " SpecMAT Triangular pad coordinates (centers) - Pad Index - (x,y)" << std::endl;
44 
45  for (index i = 0; i < 3174; ++i) {
46  coordmap << i << " ";
47  coordmap << AtPadCoord[i][3][0] << " " << AtPadCoord[i][3][1];
48  coordmap << std::endl;
49  }
50 
51  coordmap.close();
52 }
53 
54 std::vector<double> TrianglesGenerator()
55 {
56  std::vector<double> allCoordinates;
57  int nRaws = 23;
58  int firstElement = 1;
59  int elementStep = 2;
60  int nElementsInSegment = nRaws * ((2 * firstElement + (nRaws - 1) * elementStep) / 2);
61  int nSegments = 6;
62  int step = 8;
63  int counter = 0;
64  int a = 0;
65  const int arrayLength = nElementsInSegment * nSegments;
66  std::vector<double> x0(arrayLength), x1(arrayLength), x2(arrayLength), xCentre(arrayLength), x0next(arrayLength),
67  x1next(arrayLength), x2next(arrayLength);
68  std::vector<double> y0(arrayLength), y1(arrayLength), y2(arrayLength), yCentre(arrayLength), y0next(arrayLength),
69  y1next(arrayLength), y2next(arrayLength);
70  double h = 4.725;
71 
72  for (int i = 1; i < nRaws + 1; i++) {
73  a = 1 + (i - 1) * 2;
74  for (int j = 1; j < a + 1; j++) {
75  if (i == 1) {
76  x0[counter] = 0.;
77  y0[counter] = 0.;
78 
79  x1[counter] = -h / sqrt(3);
80  y1[counter] = h;
81 
82  x2[counter] = h / sqrt(3);
83  y2[counter] = h;
84 
85  xCentre[counter] = 0.;
86  yCentre[counter] = h / sqrt(3);
87  } else {
88  if (j % 2 == 0) {
89  x0[counter] = x2[counter - 1];
90  y0[counter] = y2[counter - 1];
91 
92  x1[counter] = x2[counter - a + 1];
93  y1[counter] = y2[counter - a + 1];
94 
95  x2[counter] = x1[counter - a + 1];
96  y2[counter] = y1[counter - a + 1];
97 
98  xCentre[counter] = x0[counter];
99  yCentre[counter] = y0[counter] - h / sqrt(3);
100  } else {
101  x0[counter] = x0next[i - 2] + (j - 1) * (h / sqrt(3));
102  y0[counter] = y0next[i - 2];
103 
104  x1[counter] = x1next[i - 2] + (j - 1) * (h / sqrt(3));
105  y1[counter] = y1next[i - 2];
106 
107  x2[counter] = x2next[i - 2] + (j - 1) * (h / sqrt(3));
108  y2[counter] = y2next[i - 2];
109 
110  xCentre[counter] = x0[counter];
111  yCentre[counter] = y0[counter] + h / sqrt(3);
112  }
113  }
114  counter += 1;
115  }
116 
117  x0next[i - 1] = x0[counter - a] - h / sqrt(3);
118  y0next[i - 1] = y0[counter - a] + h;
119 
120  x1next[i - 1] = x1[counter - a] - h / sqrt(3);
121  y1next[i - 1] = y1[counter - a] + h;
122 
123  x2next[i - 1] = x2[counter - a] - h / sqrt(3);
124  y2next[i - 1] = y2[counter - a] + h;
125  }
126 
127  for (int i = 1; i < nSegments; i++) {
128  for (int j = 0; j < nElementsInSegment; j++) {
129  x0[i * nElementsInSegment + j] =
130  0.5 * (x0[(i - 1) * nElementsInSegment + j] - 0) - sqrt(3) * 0.5 * (y0[(i - 1) * nElementsInSegment + j]);
131  y0[i * nElementsInSegment + j] =
132  sqrt(3) * 0.5 * (x0[(i - 1) * nElementsInSegment + j] - 0) + 0.5 * (y0[(i - 1) * nElementsInSegment + j]);
133 
134  x1[i * nElementsInSegment + j] =
135  0.5 * (x1[(i - 1) * nElementsInSegment + j] - 0) - sqrt(3) * 0.5 * (y1[(i - 1) * nElementsInSegment + j]);
136  y1[i * nElementsInSegment + j] =
137  sqrt(3) * 0.5 * (x1[(i - 1) * nElementsInSegment + j] - 0) + 0.5 * (y1[(i - 1) * nElementsInSegment + j]);
138 
139  x2[i * nElementsInSegment + j] =
140  0.5 * (x2[(i - 1) * nElementsInSegment + j] - 0) - sqrt(3) * 0.5 * (y2[(i - 1) * nElementsInSegment + j]);
141  y2[i * nElementsInSegment + j] =
142  sqrt(3) * 0.5 * (x2[(i - 1) * nElementsInSegment + j] - 0) + 0.5 * (y2[(i - 1) * nElementsInSegment + j]);
143 
144  xCentre[i * nElementsInSegment + j] = 0.5 * (xCentre[(i - 1) * nElementsInSegment + j] - 0) -
145  sqrt(3) * 0.5 * (yCentre[(i - 1) * nElementsInSegment + j]);
146  yCentre[i * nElementsInSegment + j] = sqrt(3) * 0.5 * (xCentre[(i - 1) * nElementsInSegment + j] - 0) +
147  0.5 * (yCentre[(i - 1) * nElementsInSegment + j]);
148  }
149  }
150 
151  for (int i = 0; i < arrayLength; i++) {
152  allCoordinates.push_back(x0[i]);
153  allCoordinates.push_back(y0[i]);
154  allCoordinates.push_back(x1[i]);
155  allCoordinates.push_back(y1[i]);
156  allCoordinates.push_back(x2[i]);
157  allCoordinates.push_back(y2[i]);
158  allCoordinates.push_back(xCentre[i]);
159  allCoordinates.push_back(yCentre[i]);
160  }
161  return allCoordinates;
162 } // End TrianglesGenerator()
163 
164 void AtSpecMATMap::SpecMATPadPlane()
165 {
166 
167  std::vector<double> arrayAllCoordinates2;
168  arrayAllCoordinates2 = TrianglesGenerator();
169 
170  int nElementsInSegment3;
171  nElementsInSegment3 = arrayAllCoordinates2.size() / (6 * 8);
172 
173  TH2Poly *XY;
174  XY = GetPadPlane();
175  // XY = SpecMATHistoGenerator();
176  XY->SetTitle("PadPlane X vs Y");
177  XY->SetName("XY");
178  XY->GetXaxis()->SetTitle("X (mm)");
179  XY->GetYaxis()->SetTitle("Y (mm)");
180 
181  int currentBin = 0;
182  double InitialPoInt_X = 50;
183  double InitialPoInt_Y = 50;
184  double particleEnergy = 50;
185 
186  double binXcentroid = 0;
187  double binYcentroid = 0;
188 
189  // Histogram filling example
190  XY->Fill(InitialPoInt_X, InitialPoInt_Y, particleEnergy);
191  XY->Draw("COLZ");
192 
193  // Obtain coordinates of the filled pad
194  currentBin = XY->FindBin(InitialPoInt_X, InitialPoInt_Y, 0);
195 
196  if (currentBin < nElementsInSegment3) {
197  binXcentroid = arrayAllCoordinates2[8 * (currentBin) + 6];
198  binYcentroid = arrayAllCoordinates2[8 * (currentBin) + 7];
199  } else if (currentBin >= nElementsInSegment3 && currentBin < nElementsInSegment3 * 2 - 1) {
200  binXcentroid = arrayAllCoordinates2[8 * (currentBin + 1) + 6];
201  binYcentroid = arrayAllCoordinates2[8 * (currentBin + 1) + 7];
202  } else if (currentBin >= nElementsInSegment3 * 2 - 1 && currentBin < nElementsInSegment3 * 3 - 2) {
203  binXcentroid = arrayAllCoordinates2[8 * (currentBin + 2) + 6];
204  binYcentroid = arrayAllCoordinates2[8 * (currentBin + 2) + 7];
205  } else if (currentBin >= nElementsInSegment3 * 3 - 2 && currentBin < nElementsInSegment3 * 4 - 3) {
206  binXcentroid = arrayAllCoordinates2[8 * (currentBin + 3) + 6];
207  binYcentroid = arrayAllCoordinates2[8 * (currentBin + 3) + 7];
208  } else if (currentBin >= nElementsInSegment3 * 4 - 3 && currentBin < nElementsInSegment3 * 5 - 4) {
209  binXcentroid = arrayAllCoordinates2[8 * (currentBin + 4) + 6];
210  binYcentroid = arrayAllCoordinates2[8 * (currentBin + 4) + 7];
211  } else if (currentBin >= nElementsInSegment3 * 5 - 4 && currentBin < nElementsInSegment3 * 6 - 5) {
212  binXcentroid = arrayAllCoordinates2[8 * (currentBin + 5) + 6];
213  binYcentroid = arrayAllCoordinates2[8 * (currentBin + 5) + 7];
214  }
215 
216  std::cout << "Input" << std::endl;
217  std::cout << "X: " << InitialPoInt_X << " Y: " << InitialPoInt_Y << std::endl;
218  std::cout << "\nOutput" << std::endl;
219  std::cout << "Bin: " << currentBin << " binCentroidX: " << binXcentroid << " binCentroidY: " << binYcentroid
220  << std::endl;
221 } // SpecMATPadPlane
222 
223 void AtSpecMATMap::GeneratePadPlane()
224 {
225  if (fPadPlane) {
226  LOG(error) << "Skipping generation of pad plane, it is already parsed!";
227  return;
228  }
229 
230  std::cout << " SpecMAT Map : Generating the map geometry of SpecMAT " << std::endl;
231 
232  std::vector<double> arrayAllCoord;
233  arrayAllCoord = TrianglesGenerator();
234 
235  for (int i = 0; i < fNumberPads; i++) {
236  AtPadCoord[i][0][0] = arrayAllCoord[8 * i]; // x coordinate of first vertex pad
237  AtPadCoord[i][0][1] = arrayAllCoord[8 * i + 1]; // y coordinate of first vertex pad
238  AtPadCoord[i][1][0] = arrayAllCoord[8 * i + 2]; // x coordinate of second vertex pad
239  AtPadCoord[i][1][1] = arrayAllCoord[8 * i + 3]; // y coordinate of second vertex pad
240  AtPadCoord[i][2][0] = arrayAllCoord[8 * i + 4]; // x coordinate of third vertex pad
241  AtPadCoord[i][2][1] = arrayAllCoord[8 * i + 5]; // y coordinate of third vertex pad
242  AtPadCoord[i][3][0] = arrayAllCoord[8 * i + 6]; // x coordinate of pad center
243  AtPadCoord[i][3][1] = arrayAllCoord[8 * i + 7]; // y coordinate of pad center
244  }
245 
246  std::cout << " A total of " << fNumberPads << " pads were generated " << std::endl;
247  kIsParsed = true;
248 
249  fPadPlane = new TH2Poly(); // NOLINT
250  std::vector<double> arrayAllCoordinates = TrianglesGenerator();
251  double x[3];
252  double y[3];
253  int arrayAllCoordinatesSize = arrayAllCoordinates.size();
254  for (int i2 = 0; i2 < arrayAllCoordinatesSize / 8; i2++) {
255  if (i2 != 0 && i2 != arrayAllCoordinatesSize * 1 / (6 * 8) && i2 != arrayAllCoordinatesSize * 2 / (6 * 8) &&
256  i2 != arrayAllCoordinatesSize * 3 / (6 * 8) && i2 != arrayAllCoordinatesSize * 4 / (6 * 8) &&
257  i2 != arrayAllCoordinatesSize * 5 / (6 * 8)) {
258  x[0] = arrayAllCoordinates[8 * i2];
259  y[0] = arrayAllCoordinates[8 * i2 + 1];
260  x[1] = arrayAllCoordinates[8 * i2 + 2];
261  y[1] = arrayAllCoordinates[8 * i2 + 3];
262  x[2] = arrayAllCoordinates[8 * i2 + 4];
263  y[2] = arrayAllCoordinates[8 * i2 + 5];
264  fPadPlane->AddBin(3, x, y);
265  }
266  }
267 }
268 
269 XYPoint AtSpecMATMap::CalcPadCenter(Int_t PadRef)
270 {
271  if (!kIsParsed) {
272  LOG(error) << " AtSpecMATMap::CalcPadCenter Error : Pad plane has not been generated or parsed";
273  return {-9999, -9999};
274  }
275 
276  if (PadRef == -1) { // Boost multi_array crashes with a negative index
277  LOG(debug) << " AtSpecMATMap::CalcPadCenter Error : Pad not found";
278  return {-9999, -9999};
279  }
280 
281  double x = -9999;
282  double y = -9999;
283 
284  if (PadRef < 529) {
285  x = AtPadCoord[PadRef][3][0];
286  y = AtPadCoord[PadRef][3][1];
287  } else if (PadRef >= 529 && PadRef < 529 * 2 - 1) {
288  x = AtPadCoord[PadRef + 1][3][0];
289  y = AtPadCoord[PadRef + 1][3][1];
290  } else if (PadRef >= 529 * 2 - 1 && PadRef < 529 * 3 - 2) {
291  x = AtPadCoord[PadRef + 2][3][0];
292  y = AtPadCoord[PadRef + 2][3][1];
293  } else if (PadRef >= 529 * 3 - 2 && PadRef < 529 * 4 - 3) {
294  x = AtPadCoord[PadRef + 3][3][0];
295  y = AtPadCoord[PadRef + 3][3][1];
296  } else if (PadRef >= 529 * 4 - 3 && PadRef < 529 * 5 - 4) {
297  x = AtPadCoord[PadRef + 4][3][0];
298  y = AtPadCoord[PadRef + 4][3][1];
299  } else if (PadRef >= 529 * 5 - 4 && PadRef < 529 * 6 - 5) {
300  x = AtPadCoord[PadRef + 5][3][0];
301  y = AtPadCoord[PadRef + 5][3][1];
302  }
303  // auto x = AtPadCoord[PadRef][3][0];
304  // auto y = AtPadCoord[PadRef][3][1];
305  return {x, y};
306 }
307 
308 ClassImp(AtSpecMATMap);
AtMap
Definition: AtMap.h:33
AtMap::GetPadPlane
TH2Poly * GetPadPlane()
Definition: AtMap.cxx:70
AtSpecMATMap::GeneratePadPlane
void GeneratePadPlane() override
Definition: AtSpecMATMap.cxx:223
AtSpecMATMap::AtSpecMATMap
AtSpecMATMap(Int_t fNumPads=3174)
Definition: AtSpecMATMap.cxx:27
XYPoint
ROOT::Math::XYPoint XYPoint
Definition: AtPatternCircle2D.cxx:16
AtSpecMATMap::~AtSpecMATMap
~AtSpecMATMap()
cNORMAL
constexpr auto cNORMAL
Definition: AtSpecMATMap.cxx:23
ClassImp
ClassImp(AtSpecMATMap)
AtMap::AtPadCoord
multiarray AtPadCoord
Definition: AtMap.h:41
AtMap::fNumberPads
UInt_t fNumberPads
Definition: AtMap.h:48
AtMap::kIsParsed
Bool_t kIsParsed
Definition: AtMap.h:43
y
const double * y
Definition: lmcurve.cxx:20
cYELLOW
constexpr auto cYELLOW
Definition: AtSpecMATMap.cxx:22
AtMap::index
multiarray::index index
Definition: AtMap.h:39
AtSpecMATMap::Dump
void Dump() override
Definition: AtSpecMATMap.cxx:38
AtSpecMATMap
Definition: AtSpecMATMap.h:20
TrianglesGenerator
std::vector< double > TrianglesGenerator()
Definition: AtSpecMATMap.cxx:54
AtSpecMATMap::CalcPadCenter
ROOT::Math::XYPoint CalcPadCenter(Int_t PadRef) override
Definition: AtSpecMATMap.cxx:269
AtSpecMATMap.h
XYPoint
ROOT::Math::XYPoint XYPoint
Definition: AtSpecMATMap.cxx:25
cRED
constexpr auto cRED
Definition: AtSpecMATMap.cxx:21
AtSpecMATMap::SpecMATPadPlane
void SpecMATPadPlane()
Definition: AtSpecMATMap.cxx:164
AtMap::fPadPlane
TH2Poly * fPadPlane
Definition: AtMap.h:47
Generated on Sat Nov 11 2023 15:30:38 for ATTPCROOT by   doxygen 1.8.18