-
Notifications
You must be signed in to change notification settings - Fork 102
/
R3BCalifaGeometry.cxx
409 lines (359 loc) · 14 KB
/
R3BCalifaGeometry.cxx
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
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
/******************************************************************************
* Copyright (C) 2019 GSI Helmholtzzentrum für Schwerionenforschung GmbH *
* Copyright (C) 2019-2024 Members of R3B Collaboration *
* *
* This software is distributed under the terms of the *
* GNU General Public Licence (GPL) version 3, *
* copied verbatim in the file "LICENSE". *
* *
* In applying this license GSI 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 <R3BCalifaGeometry.h>
#include <R3BLogger.h>
#include <FairLogger.h>
#include <TFile.h>
#include <TGeoManager.h>
#include <TGeoNavigator.h>
#include <TGeoVolume.h>
#include <TMath.h>
#include <TSystem.h>
#include <vector>
#include <boost/regex.hpp>
R3BCalifaGeometry* R3BCalifaGeometry::Instance()
{
// Returns singleton instance
static thread_local R3BCalifaGeometry instance;
return &instance;
}
R3BCalifaGeometry::R3BCalifaGeometry()
: TObject()
{
}
bool R3BCalifaGeometry::Init(int version)
{
if (!IsInitialize)
{
IsInitialize = kTRUE;
}
else
{
return kTRUE;
}
R3BLOG(info, "");
TString geoPath = gSystem->Getenv("VMCWORKDIR");
geoPath += "/geometry/";
switch (version)
{
case 2019:
// BARREL: demonstrator from 2019
geoPath += "califa_s444.geo.root";
fNumCrystals = 10000;
fGeometryVersion = version;
break;
case 2020:
// Half BARREL+ 6 IPHOS sectors
geoPath += "califa_v2019.11.geo.root";
fNumCrystals = 4864;
fGeometryVersion = version;
break;
case 2021:
// s455, S515, S509, S522 experiments: Half Barrel + Full IPHOS
geoPath += "califa_v2021.3.geo.root";
fNumCrystals = 4864;
fGeometryVersion = version;
break;
case 2024:
// S118, S091 experiments: Half Barrel (extended) + Full IPHOS + Full CEPA
geoPath += "califa_v2024.1.geo.root";
fNumCrystals = 5088;
fGeometryVersion = version;
break;
default:
// Full Barrel + Full IPHOS
geoPath += "califa_full.geo.root";
R3BLOG(warn,
"Unsupported geometry version: " << version << ", so standard full configuration will be used.");
// return kFALSE;
}
if (gGeoManager && strcmp(gGeoManager->GetTopVolume()->GetName(), "cave") == 0)
{
// Already set up (MC mode)
R3BLOG(info, "Using existing geometry");
// LOG(info) << "R3BCalifaGeometry::Opened geometry file " << geoPath;
fIsSimulation = kTRUE;
return kTRUE;
}
// Stand alone mode
R3BLOG(info, "Open geometry file " << geoPath << " for analysis.");
f = new TFile(geoPath, "READ");
TGeoVolume* v = dynamic_cast<TGeoVolume*>(f->Get("TOP"));
if (!v)
{
R3BLOG(error, "Could not open geometry file, No TOP volume");
return kFALSE;
}
v->SetName("cave");
if (!gGeoManager)
{
gGeoManager = new TGeoManager();
}
gGeoManager->SetTopVolume(v);
fIsSimulation = kFALSE;
return kTRUE;
}
R3BCalifaGeometry::~R3BCalifaGeometry()
{
R3BLOG(debug1, "");
if (gGeoManager)
{
delete gGeoManager;
}
if (f)
{
f->Close();
}
}
const TVector3& R3BCalifaGeometry::GetAngles(int iD)
{
static std::map<int, TVector3> cache;
Double_t local[3] = { 0, 0, 0 };
Double_t master[3];
const static TVector3 invalid(NAN, NAN, NAN);
std::string nameVolume;
if (cache.count(iD))
{
return cache[iD];
}
// SOLUTION FOR DOUBLE READING CHANNELS
if (iD > fNumCrystals / 2 && iD <= fNumCrystals)
iD = iD - fNumCrystals / 2; // for double reading crystals (crystals from 1 to max.[2544])
if (iD >= 1 && iD <= fNumCrystals / 2)
{
nameVolume = GetCrystalVolumePath(iD);
gGeoManager->CdTop();
if (gGeoManager->CheckPath(nameVolume.c_str()))
gGeoManager->cd(nameVolume.c_str());
else
{
R3BLOG(error, "Invalid crystal path: " << nameVolume);
return invalid;
}
gGeoManager->LocalToMaster(local, master);
}
else
{
R3BLOG(error, "Invalid crystalId: " << iD);
return invalid;
}
cache[iD] = master - fRefPoint;
return cache[iD];
}
void R3BCalifaGeometry::GetAngles(int iD, double* polar, double* azimuthal, double* rho)
{
auto& masterV = this->GetAngles(iD);
*polar = masterV.Theta();
*azimuthal = masterV.Phi();
*rho = masterV.Mag();
if (std::isnan(*polar) || std::isnan(*azimuthal) || std::isnan(*rho))
{
R3BLOG(error, " returns NaN");
}
}
std::string R3BCalifaGeometry::GetCrystalVolumePath(int iD)
{
int alveolusCopy = -1;
int alvType = -1;
int cryType = -1;
// SOLUTION FOR DOUBLE READING CHANNELS
if (iD > fNumCrystals / 2 && iD <= fNumCrystals)
{
iD = iD - fNumCrystals / 2; // for double reading crystals (crystals from 1 to max.[2544])
}
std::ostringstream oss;
static constexpr int kMaxNameAlv = 23;
std::string name_Alv[kMaxNameAlv] = { "01", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12",
"13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23" };
if (iD >= 1 && iD <= fNumCrystals / 2)
{
if (iD < 33)
{ // First ring in BARREL (single crystal per alveolus)
alvType = 1; // Alveolus type 1
alveolusCopy = iD - 1; // copy from 0 to 31
cryType = 1; // Only one crystal type (1)
}
else if (iD < 2337)
{ // All 4-crystals alveoli in BARREL and iPhos
alvType = (Int_t)((iD - 33) / 128) + 2; // Alveolus type (2, ..., 19)
alveolusCopy = (Int_t)((iD - 33 - (alvType - 2) * 128) / 4); // copy from 0 to 31
cryType = iD - 33 - (alvType - 2) * 128 - alveolusCopy * 4 + 1; // Four crystal types (1,2,3,4)
}
else if (iD < 2433)
{ // 3-crystals alveoli in iPhos
alvType = (Int_t)((iD - 2337) / 24) + 20; // Alveolus type (20, 21, 22, 23)
alveolusCopy = (Int_t)((iD - 2337 - (alvType - 20) * 24) / 3); // copy from 0 to 7
cryType = iD - 2337 - (alvType - 20) * 24 - alveolusCopy * 3 + 1; // Three crystal types (1,2,3)
}
else if (iD < 2465) // CEPA_CsI, four crystals alveoli (iD from 2433 to 2464)
{
alvType = 1;
alveolusCopy = (Int_t)((iD - 2433) / 4); // copy from 0 to 7
cryType = iD - 2433 - alveolusCopy * 4 + 1; // Four crystal types (1,2,3,4)
}
else if (iD < 2513) // CEPA_CsI, six crystals alveoli (iD from 2465 to 2512)
{
alvType = 2;
alveolusCopy = (Int_t)((iD - 2465) / 6); // copy from 0 to 7
cryType = iD - 2465 - alveolusCopy * 6 + 1; // Four crystal types (1,2,3,4,5,6)
}
else if (iD < 2545) // CEPA_CsI, four crystals alveoli (iD from 2513 to 2544)
{
alvType = 3;
alveolusCopy = (Int_t)((iD - 2513) / 4); // copy from 0 to 7
cryType = iD - 2513 - alveolusCopy * 4 + 1; // Four crystal types (1,2,3,4)
}
if (iD < 2433)
{
oss << "/cave_1/CalifaWorld_0/Alveolus_" << name_Alv[alvType - 1] << "_" << alveolusCopy << "/InnerAlv_"
<< name_Alv[alvType - 1] << "_0/WrapCry_" << name_Alv[alvType - 1] << "_" << cryType << "_0/Crystal_"
<< name_Alv[alvType - 1] << "_" << cryType << "_0";
}
else
{
oss << "/cave_1/CalifaWorld_0/Alveolus_CCSI_" << name_Alv[alvType - 1] << "_" << alveolusCopy
<< "/InnerAlv_CCSI_" << name_Alv[alvType - 1] << "_0/WrapCry_CCSI_" << name_Alv[alvType - 1] << "_"
<< cryType << "_0/Crystal_CCSI_" << name_Alv[alvType - 1] << "_" << cryType << "_0";
}
}
else
{
R3BLOG(error, "Invalid crystalId: " << iD);
}
return oss.str();
}
double R3BCalifaGeometry::GetDistanceThroughCrystals(TVector3& startVertex,
TVector3& direction,
TVector3* hitPos,
int* numCrystals,
int* crystalIds)
{
int maxNumCrystals = 0;
if (numCrystals != NULL && crystalIds != NULL)
{
maxNumCrystals = *numCrystals;
*numCrystals = 0;
}
TGeoNode* n;
gGeoManager->InitTrack(startVertex.X(),
startVertex.Y(),
startVertex.Z(),
direction.X() / direction.Mag(),
direction.Y() / direction.Mag(),
direction.Z() / direction.Mag());
Double_t distance = 0.;
const Double_t* pos;
bool inCrystal = false, wasInCrystal = false;
TString nodeName;
// find the distance to the next boundary and then extrapolate the current point/direction
// with this distance making sure that the boundary was crossed. Finally the goal would be to find the next state
// after crossing the boundary
while ((n = gGeoManager->FindNextBoundaryAndStep()))
{
nodeName = n->GetName();
if (inCrystal)
distance += gGeoManager->GetStep();
inCrystal = nodeName.BeginsWith("WrapCry_"); // can't enter to Crystal volume
if (inCrystal && maxNumCrystals != 0)
{
Int_t cid = GetCrystalId(gGeoManager->GetPath());
if (cid != -1 && (*numCrystals == 0 || cid != crystalIds[(*numCrystals) - 1]))
{
crystalIds[(*numCrystals)++] = cid;
maxNumCrystals--;
}
}
if (hitPos != NULL && inCrystal && !wasInCrystal)
{
pos = gGeoManager->GetCurrentPoint();
hitPos->SetXYZ(pos[0], pos[1], pos[2]);
wasInCrystal = true;
}
}
return distance;
}
// ROOT (or pyroot) is unable to automatically load dependencies of so's for whatever reason.
// so we try to force loading the lib during var initialisation.
// yes, this is ugly.
// static bool _dummy_load_libboost_regex = []() {
// char buf[100];
// snprintf(buf, 100, "%s/lib/libboost_regex.so", getenv("SIMPATH"));
// return gSystem->Load(buf);
//}();
int R3BCalifaGeometry::GetCrystalId(const std::string volumePath)
{ /*
* Yes, regex matching in the stepping action is probably slow.
* But if we cared about speed, we would not use string processing here
* in the first place.
* At least, this way we have well defined behaviour if something fails,
* instead of possible out of bounds access consequences (nasal demons etc).
* Also, it hopefully works with the 2019/s444 geo file too, now.
*/
static auto restr = "(Alveolus|Alveolus_CCSI)_([0-9]+)_([0-9]+).*(Crystal|Crystal_CCSI)_[^_]+_([0-9]+)_";
static auto re = boost::regex(restr, boost::regex::extended);
boost::cmatch m;
if (!boost::regex_search(volumePath.c_str(), m, re))
{
R3BLOG(error,
"\"" << volumePath
<< "\"\n"
"does not match RE \""
<< restr << "\".\n");
return 0;
}
int crystalId = 0;
bool isCCSI = m[1].str() == "Alveolus_CCSI"; // Adding CEPA CSI
int alvType = std::stoi(m[2].str()); // converting to int the alveolus type
int alveolusCopy = std::stoi(m[3].str()); // converting to int the alveolus copy
int cryType = std::stoi(m[5].str()); // converting to int the crystal type
bool invalid = kFALSE;
// cryType runs from 1 to 4 while alvType runs from 1 to 23, otherwise invaled
if (cryType < 1 || cryType > 4 || alvType < 1 || alvType > 23)
invalid = kTRUE;
// EXCEPTION: cryType runs up to 6 in CEPA CSI alveoli 2
if (isCCSI && alvType == 2 && (cryType == 5 || cryType == 6))
invalid = kFALSE;
if (invalid)
{
R3BLOG(error, "Wrong crystal numbers (1)");
LOG(info) << "---- cryType: " << cryType << " alvType: " << alvType;
LOG(info) << "path=" << volumePath;
return 0;
}
if (isCCSI)
{
if (alvType == 1)
crystalId = 2433 + alveolusCopy * 4 + (cryType - 1); // four crystal per alveolus (2433 to 2464)
if (alvType == 2)
crystalId = 2465 + alveolusCopy * 6 + (cryType - 1); // six crystal per alveolus (2465 to 2512)
if (alvType == 3)
crystalId = 2513 + alveolusCopy * 4 + (cryType - 1); // four crystal per alveolus (2513 to 2544)
}
else
{
if (alvType == 1)
crystalId = 1 + alveolusCopy; // first alveoli ring, one crystal per alveolus
else if (alvType < 20)
crystalId = 33 + (alvType - 2) * 128 + alveolusCopy * 4 + (cryType - 1); // four crystal per alveolus
else
crystalId = 2337 + (alvType - 20) * 24 + alveolusCopy * 3 + (cryType - 1); // three crystal per alveolus
}
if (crystalId < 1 || crystalId > fNumCrystals / 2)
{ // crystalId runs from 1 to fNumCrystals/2
R3BLOG(error, "Wrong crystal numbers (2)");
LOG(info) << "---- crystalId: " << crystalId;
return 0;
}
return crystalId;
}
ClassImp(R3BCalifaGeometry)