Project
Loading...
Searching...
No Matches
Geometry.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
16
17#include "TMath.h"
18#include <Math/Point2D.h>
19#include <Math/Vector2D.h>
20#include <ECalBase/Geometry.h>
23using namespace o2::ecal;
26
27//==============================================================================
28Geometry::Geometry()
29{
30 auto& pars = ECalBaseParam::Instance();
31 pars.updateFromFile("o2sim_configuration.ini", "ECalBase");
32 pars.printKeyValues(false, true, true, false);
33 mRMin = pars.rMin;
34 mNSuperModules = pars.nSuperModules;
35 mNCrystalModulesZ = pars.nCrystalModulesZ;
36 mNSamplingModulesZ = pars.nSamplingModulesZ;
37 mNCrystalModulesPhi = pars.nCrystalModulesPhi;
38 mNSamplingModulesPhi = pars.nSamplingModulesPhi;
39 mCrystalModW = pars.crystalModuleWidth;
40 mSamplingModW = pars.samplingModuleWidth;
41 mMarginCrystalToSampling = pars.marginCrystalToSampling;
42 mCrystalAlpha = pars.crystalAlphaDeg * TMath::DegToRad();
43 mSamplingAlpha = pars.samplingAlphaDeg * TMath::DegToRad();
44 mNModulesZ = 2 * mNSamplingModulesZ + 2 * mNCrystalModulesZ;
46}
47
48//==============================================================================
50{
51 return mNModulesZ;
52}
53
54//==============================================================================
56{
57 return mNSuperModules * (mNCrystalModulesPhi > mNSamplingModulesPhi ? mNCrystalModulesPhi : mNSamplingModulesPhi);
58}
59
60//==============================================================================
62{
63 double superModuleDeltaPhi = TwoPI / mNSuperModules;
64 double crystalDeltaPhi = getCrystalDeltaPhi();
65 return (superModuleDeltaPhi - crystalDeltaPhi * mNCrystalModulesPhi) / 2.;
66}
67
68//==============================================================================
70{
71 double superModuleDeltaPhi = TwoPI / mNSuperModules;
72 double samplingDeltaPhi = getSamplingDeltaPhi();
73 return (superModuleDeltaPhi - samplingDeltaPhi * mNSamplingModulesPhi) / 2.;
74}
75
77{
78 double theta = std::atan(mRMin / getFrontFaceZatMinR(i));
79 return -std::log(std::tan(theta / 2.));
80}
81
82//==============================================================================
84{
85 if (mFrontFaceCenterR.size() > 0) {
86 return;
87 }
88 mFrontFaceCenterTheta.resize(mNCrystalModulesZ + mNSamplingModulesZ);
89 mFrontFaceZatMinR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
90 mFrontFaceCenterR.resize(mNCrystalModulesZ + mNSamplingModulesZ);
91 mFrontFaceCenterZ.resize(mNCrystalModulesZ + mNSamplingModulesZ);
92 mTanBeta.resize(mNCrystalModulesZ + mNSamplingModulesZ);
93 mFrontFaceCenterSamplingPhi.resize(mNSuperModules * mNSamplingModulesPhi);
94 mFrontFaceCenterCrystalPhi.resize(mNSuperModules * mNCrystalModulesPhi);
95
96 double superModuleDeltaPhi = TwoPI / mNSuperModules;
97 double crystalDeltaPhi = getCrystalDeltaPhi();
98 double samplingDeltaPhi = getSamplingDeltaPhi();
99 double crystalPhiMin = getCrystalPhiMin();
100 double samplingPhiMin = getSamplingPhiMin();
101 for (int ism = 0; ism < mNSuperModules; ism++) {
102 // crystal
103 for (int i = 0; i < mNCrystalModulesPhi; i++) {
104 double phi0 = superModuleDeltaPhi * ism + crystalPhiMin + crystalDeltaPhi * i;
105 mFrontFaceCenterCrystalPhi[ism * mNCrystalModulesPhi + i] = phi0;
106 }
107 // sampling
108 for (int i = 0; i < mNSamplingModulesPhi; i++) {
109 double phi0 = superModuleDeltaPhi * ism + samplingPhiMin + samplingDeltaPhi * i;
110 mFrontFaceCenterSamplingPhi[ism * mNSamplingModulesPhi + i] = phi0;
111 }
112 }
113
114 double theta0 = PIHalf - mCrystalAlpha;
115 double zAtMinR = mCrystalModW * std::cos(mCrystalAlpha);
116
117 for (int m = 0; m < mNCrystalModulesZ; m++) {
118 mTanBeta[m] = std::sin(theta0 - mCrystalAlpha) * mCrystalModW / 2 / mRMin;
119 ROOT::Math::Polar2DVector vMid21(mCrystalModW / 2., PIHalf + theta0);
120 ROOT::Math::XYPoint pAtMinR(zAtMinR, mRMin);
121 ROOT::Math::XYPoint pc = pAtMinR + vMid21;
122 mFrontFaceZatMinR[m] = zAtMinR;
123 mFrontFaceCenterZ[m] = pc.x();
124 mFrontFaceCenterR[m] = pc.y();
125 mFrontFaceCenterTheta[m] = theta0;
126 theta0 -= 2 * mCrystalAlpha;
127 zAtMinR += mCrystalModW * std::cos(mCrystalAlpha) / std::sin(theta0 + mCrystalAlpha);
128 }
129
130 theta0 = mFrontFaceCenterTheta[mNCrystalModulesZ - 1] - mCrystalAlpha - mSamplingAlpha;
131 zAtMinR = mFrontFaceZatMinR[mNCrystalModulesZ - 1];
132 zAtMinR += mSamplingModW * std::cos(mSamplingAlpha) / std::sin(theta0 + mSamplingAlpha);
133 zAtMinR += mMarginCrystalToSampling;
134
135 for (int m = 0; m < mNSamplingModulesZ; m++) {
136 int i = m + mNCrystalModulesZ;
137 mTanBeta[i] = std::sin(theta0 - mSamplingAlpha) * mSamplingModW / 2 / mRMin;
138 ROOT::Math::Polar2DVector vMid21(mSamplingModW / 2., PIHalf + theta0);
139 ROOT::Math::XYPoint pAtMinR(zAtMinR, mRMin);
140 ROOT::Math::XYPoint pc = pAtMinR + vMid21;
141 mFrontFaceZatMinR[i] = zAtMinR;
142 mFrontFaceCenterZ[i] = pc.x();
143 mFrontFaceCenterR[i] = pc.y();
144 mFrontFaceCenterTheta[i] = theta0;
145 theta0 -= 2 * mSamplingAlpha;
146 zAtMinR += mSamplingModW * std::cos(mSamplingAlpha) / std::sin(theta0 + mSamplingAlpha);
147 }
148}
149
150int Geometry::getCellID(int moduleId, int sectorId, bool isCrystal)
151{
152 int cellID = 0;
153 if (isCrystal) {
154 if (moduleId % 2 == 0) { // crystal at positive eta
155 cellID = sectorId * mNModulesZ + moduleId / 2 + mNSamplingModulesZ + mNCrystalModulesZ;
156 } else { // crystal at negative eta
157 cellID = sectorId * mNModulesZ - moduleId / 2 + mNSamplingModulesZ + mNCrystalModulesZ - 1;
158 }
159 } else {
160 if (sectorId % 2 == 0) { // sampling at positive eta
161 cellID = sectorId / 2 * mNModulesZ + moduleId + mNSamplingModulesZ + mNCrystalModulesZ * 2;
162 } else { // sampling at negative eta
163 cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ - 1;
164 }
165 }
166 return cellID;
167}
168
169//==============================================================================
170std::pair<int, int> Geometry::globalRowColFromIndex(int cellID) const
171{
172 int ip = cellID / mNModulesZ; // row
173 int iz = cellID % mNModulesZ; // col
174 return {ip, iz};
175}
176
177//==============================================================================
178bool Geometry::isCrystal(int cellID)
179{
180 auto [row, col] = globalRowColFromIndex(cellID);
181 return (col >= mNSamplingModulesZ && col < mNSamplingModulesZ + 2 * mNCrystalModulesZ);
182}
183
184//==============================================================================
185std::pair<int, int> Geometry::getSectorChamber(int cellId) const
186{
187 int iphi = cellId / mNModulesZ;
188 int iz = cellId % mNModulesZ;
189 return getSectorChamber(iphi, iz);
190}
191
192//==============================================================================
193std::pair<int, int> Geometry::getSectorChamber(int iphi, int iz) const
194{
195 int chamber = iz < mNSamplingModulesZ ? 0 : (iz < mNSamplingModulesZ + 2 * mNCrystalModulesZ ? 1 : 2);
196 int sector = iphi / (chamber == 1 ? mNCrystalModulesPhi : mNSamplingModulesPhi);
197 return {sector, chamber};
198}
199
200//==============================================================================
201void Geometry::detIdToRelIndex(int cellId, int& chamber, int& sector, int& iphi, int& iz) const
202{
203 // 3 chambers - sampling/crystal/sampling
204 iphi = cellId / mNModulesZ;
205 iz = cellId % mNModulesZ;
206 auto pair = getSectorChamber(iphi, iz);
207 sector = pair.first;
208 chamber = pair.second;
209}
210
211//==============================================================================
212void Geometry::detIdToGlobalPosition(int detId, double& x, double& y, double& z)
213{
214 int chamber, sector, iphi, iz;
215 detIdToRelIndex(detId, chamber, sector, iphi, iz);
216 double r = 0;
217 if (iz < mNSamplingModulesZ + mNCrystalModulesZ) {
218 z = -mFrontFaceCenterZ[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
219 r = mFrontFaceCenterR[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1];
220 } else {
221 z = mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
222 r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)];
223 }
224 double phi = chamber == 1 ? mFrontFaceCenterCrystalPhi[iphi] : mFrontFaceCenterSamplingPhi[iphi];
225 x = r * std::cos(phi);
226 y = r * std::sin(phi);
227}
228
229//==============================================================================
230int Geometry::areNeighboursVertex(int detId1, int detId2) const
231{
232 int ch1, sector1, iphi1, iz1;
233 int ch2, sector2, iphi2, iz2;
234 detIdToRelIndex(detId1, ch1, sector1, iphi1, iz1);
235 detIdToRelIndex(detId2, ch2, sector2, iphi2, iz2);
236 if (sector1 != sector2 || ch1 != ch2) {
237 return 0;
238 }
239 if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1) {
240 return 1;
241 }
242 return 0;
243}
244
245//==============================================================================
246bool Geometry::isAtTheEdge(int cellId)
247{
248 auto [row, col] = globalRowColFromIndex(cellId);
249 if (col == 0) {
250 return 1;
251 } else if (col == mNSamplingModulesZ) {
252 return 1;
253 } else if (col == mNSamplingModulesZ - 1) {
254 return 1;
255 } else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ) {
256 return 1;
257 } else if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1) {
258 return 1;
259 } else if (col == mNModulesZ - 1) {
260 return 1;
261 }
262 for (int m = 0; m <= mNSuperModules; m++) {
263 if (isCrystal(cellId)) {
264 if (row == m * mNCrystalModulesPhi) {
265 return 1;
266 } else if (row == m * mNCrystalModulesPhi - 1) {
267 return 1;
268 }
269 } else {
270 if (row == m * mNSamplingModulesPhi) {
271 return 1;
272 } else if (row == m * mNSamplingModulesPhi - 1) {
273 return 1;
274 }
275 }
276 }
277 return 0;
278}
Geometry parameters configurable via o2-sim –configKeyValues.
int32_t i
useful math constants
uint32_t col
Definition RawData.h:4
Geometry helper class.
int getNcols() const
Definition Geometry.cxx:49
std::pair< int, int > globalRowColFromIndex(int cellID) const
Definition Geometry.cxx:170
double getSamplingPhiMin()
Definition Geometry.cxx:69
bool isAtTheEdge(int cellId)
Definition Geometry.cxx:246
void detIdToGlobalPosition(int detId, double &x, double &y, double &z)
Definition Geometry.cxx:212
std::pair< int, int > getSectorChamber(int cellId) const
Definition Geometry.cxx:185
bool isCrystal(int cellID)
Definition Geometry.cxx:178
double getFrontFaceZatMinR(int i)
Definition Geometry.h:51
int getNrows() const
Definition Geometry.cxx:55
void fillFrontFaceCenterCoordinates()
Definition Geometry.cxx:83
double getSamplingDeltaPhi()
Definition Geometry.h:63
double getCrystalPhiMin()
Definition Geometry.cxx:61
void detIdToRelIndex(int cellId, int &chamber, int &sector, int &iphi, int &iz) const
Definition Geometry.cxx:201
double getCrystalDeltaPhi()
Definition Geometry.h:62
int areNeighboursVertex(int detId1, int detId2) const
Definition Geometry.cxx:230
double getFrontFaceMaxEta(int i)
Definition Geometry.cxx:76
int getCellID(int moduleId, int sectorId, bool isCrystal)
Definition Geometry.cxx:150
GLint GLenum GLint x
Definition glcorearb.h:403
const GLfloat * m
Definition glcorearb.h:4066
GLboolean r
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
constexpr float TwoPI
constexpr float PIHalf
std::vector< int > row