-
Notifications
You must be signed in to change notification settings - Fork 494
Expand file tree
/
Copy pathGeometry.cxx
More file actions
246 lines (217 loc) · 8.51 KB
/
Geometry.cxx
File metadata and controls
246 lines (217 loc) · 8.51 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
#include "PHOSBase/Geometry.h"
#include "FairLogger.h"
#include "TSystem.h"
#include "TFile.h"
using namespace o2::phos;
ClassImp(Geometry);
// module numbering:
// start from module 0 (non-existing), 1 (half-module), 2 (bottom),... 4(highest)
// absId:
// start from 1 till 4*64*56. Numbering in each module starts at bottom left and first go in z direction:
// 56 112 3584
// ... ... ...
// 1 57 ...3529
// relid[3]: (module number[1...4], iphi[1...64], iz[1...56])
//
// Then TRU channels go from 1 to 112 per branch, 2 branches per ddl
// absId = getTotalNCells() + TRUabsId ;
// relId for TRU
// relid: [DDL id=0..13] [x in 2x2 system: 1..8] [z in 2x2 system 1..28] TODO: verify with real TRU data!!!
// these initialisations are needed for a singleton
Geometry* Geometry::sGeom = nullptr;
Geometry::Geometry(const std::string_view name) : mGeoName(name)
{
std::string p = gSystem->Getenv("O2_ROOT");
p += "/share/Detectors/PHOS/files/alignment.root";
TFile fin(p.data());
//try reading rotation mathices
for (int m = 1; m < 5; m++) {
mPHOS[m] = *static_cast<TGeoHMatrix*>(fin.Get(Form("Module%d", m)));
}
fin.Close();
}
short Geometry::relToAbsId(char moduleNumber, int strip, int cell)
{
// calculates absolute cell Id from moduleNumber, strip (number) and cell (number)
// PHOS layout parameters:
const short nStrpZ = 28; // Number of strips along z-axis
const short nCrystalsInModule = 56 * 64; // Total number of crystals in module
const short nCellsXInStrip = 8; // Number of crystals in strip unit along x-axis
const short nZ = 56; // nStripZ * nCellsZInStrip
short row = nStrpZ - (strip - 1) % nStrpZ;
short col = (int)std::ceil((float)strip / (nStrpZ)) - 1;
return (moduleNumber - 1) * nCrystalsInModule + row * 2 + (col * nCellsXInStrip + (cell - 1) / 2) * nZ -
(cell & 1 ? 1 : 0);
}
bool Geometry::absToRelNumbering(short absId, char* relid)
{
// Converts the absolute numbering into the following array
// relid[0] = PHOS Module number 1:fNModules
// relid[1] = Row number inside a PHOS module (Z coordinate)
// relid[2] = Column number inside a PHOS module (Phi coordinate)
const short nZ = 56; // nStripZ * nCellsZInStrip
const short nPhi = 64; // nStripZ * nCellsZInStrip
absId--;
short phosmodulenumber = absId / (nZ * nPhi);
relid[0] = phosmodulenumber + 1;
absId -= phosmodulenumber * nPhi * nZ;
relid[1] = 1 + absId / nZ;
relid[2] = absId - (relid[1] - 1) * nZ + 1;
return true;
}
bool Geometry::truAbsToRelNumbering(short truId, char* relid)
{
//convert trigger cell Id to
truId--;
relid[0] = truId / 224; //2*112 channels // DDL id
truId = truId % 224;
relid[1] = 1 + truId % 8; // x index in TRU internal 2x2 coordinate system
relid[2] = 1 + truId / 8; // z index in TRU internal 2x2 coordinate system
return true;
}
short Geometry::truRelToAbsNumbering(const char* relId)
{
return relId[0] * 224 + // the offset of PHOS modules
relId[1] - 1 + // the offset along phi
(relId[2] - 1) * 8; // the offset along z
}
bool Geometry::truRelId2RelId(const char* truRelId, char* relId)
{
relId[0] = 1 + (truRelId[0] + 2) / 4;
relId[1] = ((truRelId[0] + 2) % 4) * 16 + truRelId[1] * 2 - 1;
relId[2] = truRelId[2] * 2 - 1;
return true;
}
short Geometry::relPosToTruId(char mod, float x, float z, short& ddl)
{
//tranform local cluster coordinates to truId
char relid[3] = {mod, static_cast<char>(ceil(x / CELLSTEP + 32.5)), static_cast<char>(ceil(z / CELLSTEP + 28.5))};
ddl = (mod - 1) * 4 + relid[1] / 16 - 2;
char truid[3] = {static_cast<char>(ddl), static_cast<char>(1 + ((relid[1] - 1) % 16) / 2), static_cast<char>(1 + (relid[2] - 1) / 2)};
return truRelToAbsNumbering(truid);
}
char Geometry::absIdToModule(short absId)
{
const short nZ = 56;
const short nPhi = 64;
return 1 + (absId - 1) / (nZ * nPhi);
}
int Geometry::areNeighbours(short absId1, short absId2)
{
// Gives the neighbourness of two digits = 0 are not neighbour but continue searching
// = 1 are neighbour
// = 2 are not neighbour but do not continue searching
// =-1 are not neighbour, continue searching, but do not look before d2 next
// time
// neighbours are defined as digits having at least a common vertex
// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
// which is compared to a digit (d2) not yet in a cluster
char relid1[3];
absToRelNumbering(absId1, relid1);
char relid2[3];
absToRelNumbering(absId2, relid2);
if (relid1[0] == relid2[0]) { // inside the same PHOS module
char rowdiff = TMath::Abs(relid1[1] - relid2[1]);
char coldiff = TMath::Abs(relid1[2] - relid2[2]);
if ((coldiff <= 1) && (rowdiff <= 1)) { // At least common vertex
return 1;
} else {
if ((relid2[1] > relid1[1]) && (relid2[2] > relid1[2] + 1)) {
return 2; // Difference in row numbers is too large to look further
}
}
return 0;
} else {
if (relid1[0] > relid2[0]) { // we switched to the next module
return -1;
}
return 2;
}
return 0;
}
void Geometry::absIdToRelPosInModule(short absId, float& x, float& z)
{
char relid[3];
absToRelNumbering(absId, relid);
x = (relid[1] - 32 - 0.5) * CELLSTEP;
z = (relid[2] - 28 - 0.5) * CELLSTEP;
}
bool Geometry::relToAbsNumbering(const char* relId, short& absId)
{
const short nZ = 56; // nStripZ * nCellsZInStrip
const short nPhi = 64; // nStripZ * nCellsZInStrip
absId =
(relId[0] - 1) * nPhi * nZ + // the offset of PHOS modules
(relId[1] - 1) * nZ + // the offset along phi
relId[2]; // the offset along z
return true;
}
//local position to absId
void Geometry::relPosToAbsId(char module, float x, float z, short& absId)
{
char relid[3] = {module, static_cast<char>(ceil(x / CELLSTEP + 32.5)), static_cast<char>(ceil(z / CELLSTEP + 28.5))};
relToAbsNumbering(relid, absId);
}
void Geometry::relPosToRelId(short module, float x, float z, char* relId)
{
relId[0] = module;
relId[1] = static_cast<char>(ceil(x / CELLSTEP + 32.5));
relId[2] = static_cast<char>(ceil(z / CELLSTEP + 28.5));
}
// convert local position in module to global position in ALICE
void Geometry::local2Global(char module, float x, float z, TVector3& globaPos) const
{
// constexpr float shiftY=-10.76; Run2
constexpr float shiftY = -1.26; //Depth-optimized
Double_t posL[3] = {x, z, shiftY};
Double_t posG[3];
mPHOS[module].LocalToMaster(posL, posG);
globaPos.SetXYZ(posG[0], posG[1], posG[2]);
}
bool Geometry::impactOnPHOS(const TVector3& vtx, const TVector3& p,
short& module, float& z, float& x) const
{
// calculates the impact coordinates on PHOS of a neutral particle
// emitted in the vertex vtx with 3-momentum p
constexpr float shiftY = -1.26; //Depth-optimized
constexpr float moduleXhalfSize = 72.16; // 18.04 / 2 * 8
constexpr float moduleZhalfSize = 64.14; // 4.51 / 2 * 28
for (short mod = 1; mod < 5; mod++) {
//create vector from (0,0,0) to center of crystal surface of imod module
double tmp[3] = {0., 0., shiftY};
double posG[3] = {0., 0., 0.};
mPHOS[mod].LocalToMaster(tmp, posG);
TVector3 n(posG[0], posG[1], posG[2]);
double direction = n.Dot(p);
if (direction <= 0.) {
continue; //momentum directed FROM module
}
double fr = (n.Mag2() - n.Dot(vtx)) / direction;
//Calculate direction in module plane
n -= vtx + fr * p;
n *= -1.;
if (TMath::Abs(n.Z()) < moduleZhalfSize && n.Pt() < moduleXhalfSize) {
module = mod;
z = n.Z();
x = TMath::Sign(n.Pt(), n.X());
//no need to return to local system since we calcilated distance from module center
//and tilts can not be significant.
return true;
}
}
//Not in acceptance
x = 0;
z = 0;
module = 0;
return false;
}