Skip to content

Commit c4f7253

Browse files
mslupeckmslupeck
authored
[FV0] Added getter and utility functions to fv0::Geometry to be used … (#3463)
1 parent 80480bd commit c4f7253

6 files changed

Lines changed: 151 additions & 16 deletions

File tree

‎Detectors/FIT/FV0/base/CMakeLists.txt‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,7 +10,9 @@
1010

1111
o2_add_library(FV0Base
1212
SOURCES src/Geometry.cxx
13-
PUBLIC_LINK_LIBRARIES ROOT::Geom FairRoot::Base)
13+
PUBLIC_LINK_LIBRARIES ROOT::Geom
14+
FairRoot::Base
15+
O2::SimulationDataFormat)
1416

1517
o2_target_root_dictionary(FV0Base HEADERS include/FV0Base/Geometry.h)
1618

‎Detectors/FIT/FV0/base/include/FV0Base/Geometry.h‎

Lines changed: 47 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -18,10 +18,12 @@
1818
#define ALICEO2_FV0_GEOMETRY_H_
1919

2020
#include <vector>
21+
#include <array>
2122

2223
#include <TGeoMatrix.h>
2324
#include <TGeoVolume.h>
2425
#include <TVirtualMC.h>
26+
#include "MathUtils/Cartesian3D.h"
2527

2628
namespace o2
2729
{
@@ -32,9 +34,13 @@ class Geometry
3234
{
3335
public:
3436
/// Geometry type options possible to be initialized. The type of the geometry will specify which components are
35-
/// created.
37+
/// created. Geometry types
38+
/// -> eUnitialized => no parts
39+
/// -> eOnlySensitive => only sensitive detector parts
40+
/// -> eRough => sensitive parts and rough structural elements
41+
/// -> eFull => complete, detailed geometry (including screws, etc.)
3642
enum EGeoType {
37-
eUninitilized,
43+
eUninitialized,
3844
eOnlySensitive,
3945
eRough,
4046
eFull
@@ -53,21 +59,15 @@ class Geometry
5359
/// Default constructor.
5460
/// It must be kept public for root persistency purposes,
5561
/// but should never be called by the outside world
56-
Geometry() : mGeometryType(eUninitilized), mLeftTransformation(nullptr), mRightTransformation(nullptr){};
57-
58-
/// Standard constructor
59-
/// \param initType[in] The type of geometry, that will be initialized
60-
/// -> initType == eUnitialized => no parts
61-
/// -> initType == eOnlySensitive => only sensitive detector parts
62-
/// -> initType == eRough => sensitive parts and rough structural elements
63-
/// -> initType == eFull => complete, detailed geometry (including screws, etc.)
64-
/// \return -
65-
explicit Geometry(EGeoType initType);
62+
Geometry() : mGeometryType(eUninitialized), mLeftTransformation(nullptr), mRightTransformation(nullptr){};
6663

6764
/// Copy constructor.
6865
Geometry(const Geometry& geometry);
6966

70-
/// Destructor
67+
/// Access to geometry instance
68+
/// \param initType The geometry type to be initialized - if the geometry already exists this parameter is ignored
69+
static Geometry* instance(EGeoType initType = eUninitialized);
70+
7171
~Geometry();
7272

7373
/// Get the unique ID of the current scintillator cell during simulation.
@@ -91,7 +91,24 @@ class Geometry
9191
/// Build the geometry.
9292
void buildGeometry() const;
9393

94+
/// Utility functions to be accessed externally
95+
96+
/// Sets the input parameters to the position of the geometrical center of sensitive detector
97+
/// \param x x [cm].
98+
/// \param y y [cm].
99+
/// \param z z [cm].
100+
void getGlobalPosition(float& x, float& y, float& z);
101+
Point3D<float>& getCellCenter(UInt_t cellId);
102+
Point3D<float>& getReadoutCenter(UInt_t cellId);
103+
104+
/// Helper function to check if the cellId belongs to ring 5.
105+
/// \param cellId Id of the cell in range from 0 to 39.
106+
/// \return True if cellId belongs to ring 5.
107+
bool isRing5(UInt_t cellId);
108+
94109
private:
110+
explicit Geometry(EGeoType initType);
111+
95112
inline static const std::string sDetectorName = "FV0";
96113

97114
// General geometry constants
@@ -110,6 +127,13 @@ class Geometry
110127
/// Cell and scintillator constants
111128
static constexpr int sNumberOfCellSectors = 4; ///< Number of cell sectors for one half of the detector
112129
static constexpr int sNumberOfCellRings = 5; ///< Number of cell rings
130+
static constexpr int sNumberOfCells = sNumberOfCellRings * sNumberOfCellSectors * 2; ///< Number of cells
131+
static constexpr int sNumberOfReadoutChannels = sNumberOfCells + sNumberOfCellSectors * 2; ///< Number of ch (2 ch per cell in r5)
132+
133+
/// Look-up tables converting cellId to the ring and sector number
134+
static constexpr int sCellToRing[sNumberOfCells] = {0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4};
135+
static constexpr int sCellToSector[sNumberOfCells] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3, 4, 5, 6, 7};
136+
113137
/// Average cell ring radii.
114138
static constexpr float sCellRingRadii[sNumberOfCellRings + 1]{4.01, 7.3, 12.9, 21.25, 38.7, 72.115};
115139
static constexpr char sCellTypes[sNumberOfCellSectors]{'a', 'b', 'b', 'a'}; ///< Ordered cell types per half a ring
@@ -386,6 +410,10 @@ class Geometry
386410
/// \return The volume name.
387411
const std::string createVolumeName(const std::string& volumeType, int number = -1) const;
388412

413+
/// Utility methods
414+
void initializeCellCenters(); ///< To be called in constructor to initialize mCellCenter
415+
void initializeReadoutCenters(); ///< To be called in constructor to initialize mReadoutCenter
416+
389417
std::vector<std::string> mSensitiveVolumeNames; ///< The names of all the sensitive volumes
390418

391419
/// Average ring radii
@@ -430,6 +458,12 @@ class Geometry
430458
TGeoMatrix* mLeftTransformation; ///< Transformation for the left part of the detector
431459
TGeoMatrix* mRightTransformation; ///< Transformation for the right part of the detector
432460

461+
/// Utility arrays derived from constants
462+
std::array<Point3D<float>, sNumberOfCells> mCellCenter; ///< Center of each scintillator cell
463+
std::array<Point3D<float>, sNumberOfReadoutChannels> mReadoutCenter; ///< Similar to mCellCenter, cells in r5 are additionally divided
464+
465+
static Geometry* sInstance; ///< \brief Singleton instance
466+
433467
ClassDefNV(Geometry, 1);
434468
};
435469
} // namespace fv0

‎Detectors/FIT/FV0/base/src/Geometry.cxx‎

Lines changed: 78 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,11 @@ using namespace o2::fv0;
3535

3636
Geometry::Geometry(EGeoType initType) : mGeometryType(initType)
3737
{
38-
initializeGeometry();
38+
initializeCellCenters();
39+
initializeReadoutCenters();
40+
if (initType != eUninitialized) {
41+
initializeGeometry();
42+
}
3943
}
4044

4145
Geometry::Geometry(const Geometry& geometry) : mGeometryType(geometry.mGeometryType), mLeftTransformation(nullptr), mRightTransformation(nullptr)
@@ -105,6 +109,28 @@ void Geometry::buildGeometry() const
105109
vALIC->AddNode(vFV0, 1, new TGeoTranslation(sXGlobal, sYGlobal, sZGlobal));
106110
}
107111

112+
void Geometry::getGlobalPosition(float& x, float& y, float& z)
113+
{
114+
x = sXGlobal;
115+
y = sYGlobal;
116+
z = sZGlobal;
117+
}
118+
119+
Point3D<float>& Geometry::getCellCenter(UInt_t cellId)
120+
{
121+
return mCellCenter.at(cellId);
122+
}
123+
124+
Point3D<float>& Geometry::getReadoutCenter(UInt_t cellId)
125+
{
126+
return mReadoutCenter.at(cellId);
127+
}
128+
129+
bool Geometry::isRing5(UInt_t cellId)
130+
{
131+
return cellId >= (sNumberOfCellRings - 1) * sNumberOfCellSectors * 2;
132+
}
133+
108134
void Geometry::initializeGeometry()
109135
{
110136
initializeMaps();
@@ -1185,3 +1211,54 @@ const std::string Geometry::createVolumeName(const std::string& volumeType, cons
11851211
{
11861212
return sDetectorName + volumeType + ((number >= 0) ? std::to_string(number) : "");
11871213
}
1214+
1215+
void Geometry::initializeCellCenters()
1216+
{
1217+
const float phi0 = 67.5 * TMath::DegToRad(); // starting phi of one of the sectors
1218+
const float dphi = 45. * TMath::DegToRad(); // phi difference between neighbouring sectors
1219+
const float lutSect2Phi[sNumberOfCellSectors * 2] = {phi0, phi0 - dphi, phi0 - 2 * dphi, phi0 - 3 * dphi, phi0 + dphi, phi0 + 2 * dphi, phi0 + 3 * dphi, phi0 + 4 * dphi};
1220+
for (int cellId = 0; cellId < sNumberOfCells; cellId++) {
1221+
float r = 0.5 * (sCellRingRadii[sCellToRing[cellId]] + sCellRingRadii[sCellToRing[cellId] + 1]);
1222+
double x = sXGlobal + r * TMath::Cos(lutSect2Phi[sCellToSector[cellId]]);
1223+
double y = sYGlobal + r * TMath::Sin(lutSect2Phi[sCellToSector[cellId]]);
1224+
1225+
Point3D<float>* p = &mCellCenter.at(cellId);
1226+
p->SetCoordinates(x, y, sZGlobal);
1227+
}
1228+
}
1229+
1230+
void Geometry::initializeReadoutCenters()
1231+
{
1232+
for (int channelId = 0; channelId < sNumberOfReadoutChannels; channelId++) {
1233+
Point3D<float>* p = &mReadoutCenter.at(channelId);
1234+
if (!isRing5(channelId)) {
1235+
p->SetCoordinates(getCellCenter(channelId).X(), getCellCenter(channelId).Y(), getCellCenter(channelId).Z());
1236+
} else {
1237+
const int numberOfSectorsR5 = sNumberOfCellSectors * 4; // from both halves of the detector
1238+
const float phi0 = 78.75 * TMath::DegToRad(); // starting phi of one of the sectors
1239+
const float dphi = 22.5 * TMath::DegToRad(); // phi difference between neighbouring sectors
1240+
const float lutReadoutSect2Phi[numberOfSectorsR5] =
1241+
{phi0 - 0 * dphi, phi0 - 1 * dphi, phi0 - 2 * dphi, phi0 - 3 * dphi,
1242+
phi0 - 4 * dphi, phi0 - 5 * dphi, phi0 - 6 * dphi, phi0 - 7 * dphi,
1243+
phi0 + 1 * dphi, phi0 + 2 * dphi, phi0 + 3 * dphi, phi0 + 4 * dphi,
1244+
phi0 + 5 * dphi, phi0 + 6 * dphi, phi0 + 7 * dphi, phi0 + 8 * dphi};
1245+
1246+
int iReadoutSector = channelId - ((sNumberOfCellRings - 1) * sNumberOfCellSectors * 2);
1247+
float r = 0.5 * (sCellRingRadii[4] + sCellRingRadii[5]);
1248+
double x = sXGlobal + r * TMath::Cos(lutReadoutSect2Phi[iReadoutSector]);
1249+
double y = sYGlobal + r * TMath::Sin(lutReadoutSect2Phi[iReadoutSector]);
1250+
p->SetCoordinates(x, y, sZGlobal);
1251+
}
1252+
}
1253+
}
1254+
1255+
Geometry* Geometry::sInstance = nullptr;
1256+
1257+
//Singleton access
1258+
Geometry* Geometry::instance(EGeoType initType)
1259+
{
1260+
if (!sInstance)
1261+
LOG(INFO) << "FV0 geometry instance created";
1262+
sInstance = new Geometry(initType);
1263+
return sInstance;
1264+
}

‎Detectors/FIT/FV0/simulation/include/FV0Simulation/Digitizer.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,7 @@ class Digitizer
9292
static Double_t PmtResponse(Double_t x);
9393
static Double_t PmtResponse(Double_t* x, Double_t*);
9494
static Double_t SinglePhESpectrum(Double_t* x, Double_t* par);
95+
float getDistFromCellCenter(UInt_t cellId, double hitx, double hity);
9596

9697
ClassDefNV(Digitizer, 1);
9798
};

‎Detectors/FIT/FV0/simulation/src/Detector.cxx‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -267,7 +267,7 @@ void Detector::ConstructGeometry()
267267
{
268268
LOG(INFO) << "FV0: Constructing geometry";
269269
createMaterials();
270-
mGeometry = new Geometry(Geometry::eFull);
270+
mGeometry = Geometry::instance(Geometry::eFull);
271271
// mGeometry->enableComponent(Geometry::eScintillator, false);
272272
// mGeometry->enableComponent(Geometry::ePlastics, false);
273273
// mGeometry->enableComponent(Geometry::eFibers, false);

‎Detectors/FIT/FV0/simulation/src/Digitizer.cxx‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99
// or submit itself to any jurisdiction.
1010

1111
#include "FV0Simulation/Digitizer.h"
12+
#include "FV0Base/Geometry.h"
1213

1314
#include <TRandom.h>
1415
#include <algorithm>
@@ -122,6 +123,9 @@ void Digitizer::process(const std::vector<o2::fv0::Hit>& hits)
122123
Float_t const t = hit.GetTime() * 1e9 + FV0DigParam::Instance().pmtTransitTime;
123124
Float_t const charge = TMath::Qe() * FV0DigParam::Instance().pmtGain * mBinSize / mPmtTimeIntegral;
124125

126+
// Example how to access the fields necessary to split the readout of cells in ring5 to two PMTs
127+
// LOG(INFO) << detId << " " << Geometry::instance()->isRing5(detId) << " " << hit.GetX() << ", " << hit.GetY() << " " << getDistFromCellCenter(detId, hit.GetX(), hit.GetY());
128+
125129
auto& analogSignal = mPmtChargeVsTime[detId];
126130

127131
for (Int_t iPhE = 0; iPhE < nPhE; ++iPhE) {
@@ -254,3 +258,20 @@ Double_t Digitizer::SinglePhESpectrum(Double_t* x, Double_t*)
254258
return (TMath::Poisson(x[0], FV0DigParam::Instance().pmtNbOfSecElec) +
255259
FV0DigParam::Instance().pmtTransparency * TMath::Poisson(x[0], 1.0));
256260
}
261+
262+
// The Distance is positive for top half-sectors (when the hit position is above the cell center (has higher y))
263+
// TODO: performance check needed
264+
float Digitizer::getDistFromCellCenter(UInt_t cellId, double hitx, double hity)
265+
{
266+
Geometry* geo = Geometry::instance();
267+
268+
// Parametrize the line (ax+by+c=0) that crosses the detector center and the cell's middle point
269+
Point3D<float>* pCell = &geo->getCellCenter(cellId);
270+
float x0, y0, z0;
271+
geo->getGlobalPosition(x0, y0, z0);
272+
double a = -(y0 - pCell->Y()) / (x0 - pCell->X());
273+
double b = 1;
274+
double c = -(y0 - a * x0);
275+
// Return the distance from hit to this line
276+
return (a * hitx + b * hity + c) / TMath::Sqrt(a * a + b * b);
277+
}

0 commit comments

Comments
 (0)