diff --git a/Detectors/FIT/FV0/base/CMakeLists.txt b/Detectors/FIT/FV0/base/CMakeLists.txt index cae145bab2700..d96cb91315d66 100644 --- a/Detectors/FIT/FV0/base/CMakeLists.txt +++ b/Detectors/FIT/FV0/base/CMakeLists.txt @@ -10,7 +10,9 @@ o2_add_library(FV0Base SOURCES src/Geometry.cxx - PUBLIC_LINK_LIBRARIES ROOT::Geom FairRoot::Base) + PUBLIC_LINK_LIBRARIES ROOT::Geom + FairRoot::Base + O2::SimulationDataFormat) o2_target_root_dictionary(FV0Base HEADERS include/FV0Base/Geometry.h) diff --git a/Detectors/FIT/FV0/base/include/FV0Base/Geometry.h b/Detectors/FIT/FV0/base/include/FV0Base/Geometry.h index 1472647e3a6bd..97e9d02cbdbc6 100644 --- a/Detectors/FIT/FV0/base/include/FV0Base/Geometry.h +++ b/Detectors/FIT/FV0/base/include/FV0Base/Geometry.h @@ -18,10 +18,12 @@ #define ALICEO2_FV0_GEOMETRY_H_ #include +#include #include #include #include +#include "MathUtils/Cartesian3D.h" namespace o2 { @@ -32,9 +34,13 @@ class Geometry { public: /// Geometry type options possible to be initialized. The type of the geometry will specify which components are - /// created. + /// created. Geometry types + /// -> eUnitialized => no parts + /// -> eOnlySensitive => only sensitive detector parts + /// -> eRough => sensitive parts and rough structural elements + /// -> eFull => complete, detailed geometry (including screws, etc.) enum EGeoType { - eUninitilized, + eUninitialized, eOnlySensitive, eRough, eFull @@ -53,21 +59,15 @@ class Geometry /// Default constructor. /// It must be kept public for root persistency purposes, /// but should never be called by the outside world - Geometry() : mGeometryType(eUninitilized), mLeftTransformation(nullptr), mRightTransformation(nullptr){}; - - /// Standard constructor - /// \param initType[in] The type of geometry, that will be initialized - /// -> initType == eUnitialized => no parts - /// -> initType == eOnlySensitive => only sensitive detector parts - /// -> initType == eRough => sensitive parts and rough structural elements - /// -> initType == eFull => complete, detailed geometry (including screws, etc.) - /// \return - - explicit Geometry(EGeoType initType); + Geometry() : mGeometryType(eUninitialized), mLeftTransformation(nullptr), mRightTransformation(nullptr){}; /// Copy constructor. Geometry(const Geometry& geometry); - /// Destructor + /// Access to geometry instance + /// \param initType The geometry type to be initialized - if the geometry already exists this parameter is ignored + static Geometry* instance(EGeoType initType = eUninitialized); + ~Geometry(); /// Get the unique ID of the current scintillator cell during simulation. @@ -91,7 +91,24 @@ class Geometry /// Build the geometry. void buildGeometry() const; + /// Utility functions to be accessed externally + + /// Sets the input parameters to the position of the geometrical center of sensitive detector + /// \param x x [cm]. + /// \param y y [cm]. + /// \param z z [cm]. + void getGlobalPosition(float& x, float& y, float& z); + Point3D& getCellCenter(UInt_t cellId); + Point3D& getReadoutCenter(UInt_t cellId); + + /// Helper function to check if the cellId belongs to ring 5. + /// \param cellId Id of the cell in range from 0 to 39. + /// \return True if cellId belongs to ring 5. + bool isRing5(UInt_t cellId); + private: + explicit Geometry(EGeoType initType); + inline static const std::string sDetectorName = "FV0"; // General geometry constants @@ -110,6 +127,13 @@ class Geometry /// Cell and scintillator constants static constexpr int sNumberOfCellSectors = 4; ///< Number of cell sectors for one half of the detector static constexpr int sNumberOfCellRings = 5; ///< Number of cell rings + static constexpr int sNumberOfCells = sNumberOfCellRings * sNumberOfCellSectors * 2; ///< Number of cells + static constexpr int sNumberOfReadoutChannels = sNumberOfCells + sNumberOfCellSectors * 2; ///< Number of ch (2 ch per cell in r5) + + /// Look-up tables converting cellId to the ring and sector number + 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}; + 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}; + /// Average cell ring radii. static constexpr float sCellRingRadii[sNumberOfCellRings + 1]{4.01, 7.3, 12.9, 21.25, 38.7, 72.115}; static constexpr char sCellTypes[sNumberOfCellSectors]{'a', 'b', 'b', 'a'}; ///< Ordered cell types per half a ring @@ -386,6 +410,10 @@ class Geometry /// \return The volume name. const std::string createVolumeName(const std::string& volumeType, int number = -1) const; + /// Utility methods + void initializeCellCenters(); ///< To be called in constructor to initialize mCellCenter + void initializeReadoutCenters(); ///< To be called in constructor to initialize mReadoutCenter + std::vector mSensitiveVolumeNames; ///< The names of all the sensitive volumes /// Average ring radii @@ -430,6 +458,12 @@ class Geometry TGeoMatrix* mLeftTransformation; ///< Transformation for the left part of the detector TGeoMatrix* mRightTransformation; ///< Transformation for the right part of the detector + /// Utility arrays derived from constants + std::array, sNumberOfCells> mCellCenter; ///< Center of each scintillator cell + std::array, sNumberOfReadoutChannels> mReadoutCenter; ///< Similar to mCellCenter, cells in r5 are additionally divided + + static Geometry* sInstance; ///< \brief Singleton instance + ClassDefNV(Geometry, 1); }; } // namespace fv0 diff --git a/Detectors/FIT/FV0/base/src/Geometry.cxx b/Detectors/FIT/FV0/base/src/Geometry.cxx index 8a1fee4fed96d..ccd966048d8d3 100644 --- a/Detectors/FIT/FV0/base/src/Geometry.cxx +++ b/Detectors/FIT/FV0/base/src/Geometry.cxx @@ -35,7 +35,11 @@ using namespace o2::fv0; Geometry::Geometry(EGeoType initType) : mGeometryType(initType) { - initializeGeometry(); + initializeCellCenters(); + initializeReadoutCenters(); + if (initType != eUninitialized) { + initializeGeometry(); + } } Geometry::Geometry(const Geometry& geometry) : mGeometryType(geometry.mGeometryType), mLeftTransformation(nullptr), mRightTransformation(nullptr) @@ -105,6 +109,28 @@ void Geometry::buildGeometry() const vALIC->AddNode(vFV0, 1, new TGeoTranslation(sXGlobal, sYGlobal, sZGlobal)); } +void Geometry::getGlobalPosition(float& x, float& y, float& z) +{ + x = sXGlobal; + y = sYGlobal; + z = sZGlobal; +} + +Point3D& Geometry::getCellCenter(UInt_t cellId) +{ + return mCellCenter.at(cellId); +} + +Point3D& Geometry::getReadoutCenter(UInt_t cellId) +{ + return mReadoutCenter.at(cellId); +} + +bool Geometry::isRing5(UInt_t cellId) +{ + return cellId >= (sNumberOfCellRings - 1) * sNumberOfCellSectors * 2; +} + void Geometry::initializeGeometry() { initializeMaps(); @@ -1185,3 +1211,54 @@ const std::string Geometry::createVolumeName(const std::string& volumeType, cons { return sDetectorName + volumeType + ((number >= 0) ? std::to_string(number) : ""); } + +void Geometry::initializeCellCenters() +{ + const float phi0 = 67.5 * TMath::DegToRad(); // starting phi of one of the sectors + const float dphi = 45. * TMath::DegToRad(); // phi difference between neighbouring sectors + 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}; + for (int cellId = 0; cellId < sNumberOfCells; cellId++) { + float r = 0.5 * (sCellRingRadii[sCellToRing[cellId]] + sCellRingRadii[sCellToRing[cellId] + 1]); + double x = sXGlobal + r * TMath::Cos(lutSect2Phi[sCellToSector[cellId]]); + double y = sYGlobal + r * TMath::Sin(lutSect2Phi[sCellToSector[cellId]]); + + Point3D* p = &mCellCenter.at(cellId); + p->SetCoordinates(x, y, sZGlobal); + } +} + +void Geometry::initializeReadoutCenters() +{ + for (int channelId = 0; channelId < sNumberOfReadoutChannels; channelId++) { + Point3D* p = &mReadoutCenter.at(channelId); + if (!isRing5(channelId)) { + p->SetCoordinates(getCellCenter(channelId).X(), getCellCenter(channelId).Y(), getCellCenter(channelId).Z()); + } else { + const int numberOfSectorsR5 = sNumberOfCellSectors * 4; // from both halves of the detector + const float phi0 = 78.75 * TMath::DegToRad(); // starting phi of one of the sectors + const float dphi = 22.5 * TMath::DegToRad(); // phi difference between neighbouring sectors + const float lutReadoutSect2Phi[numberOfSectorsR5] = + {phi0 - 0 * dphi, phi0 - 1 * dphi, phi0 - 2 * dphi, phi0 - 3 * dphi, + phi0 - 4 * dphi, phi0 - 5 * dphi, phi0 - 6 * dphi, phi0 - 7 * dphi, + phi0 + 1 * dphi, phi0 + 2 * dphi, phi0 + 3 * dphi, phi0 + 4 * dphi, + phi0 + 5 * dphi, phi0 + 6 * dphi, phi0 + 7 * dphi, phi0 + 8 * dphi}; + + int iReadoutSector = channelId - ((sNumberOfCellRings - 1) * sNumberOfCellSectors * 2); + float r = 0.5 * (sCellRingRadii[4] + sCellRingRadii[5]); + double x = sXGlobal + r * TMath::Cos(lutReadoutSect2Phi[iReadoutSector]); + double y = sYGlobal + r * TMath::Sin(lutReadoutSect2Phi[iReadoutSector]); + p->SetCoordinates(x, y, sZGlobal); + } + } +} + +Geometry* Geometry::sInstance = nullptr; + +//Singleton access +Geometry* Geometry::instance(EGeoType initType) +{ + if (!sInstance) + LOG(INFO) << "FV0 geometry instance created"; + sInstance = new Geometry(initType); + return sInstance; +} diff --git a/Detectors/FIT/FV0/simulation/include/FV0Simulation/Digitizer.h b/Detectors/FIT/FV0/simulation/include/FV0Simulation/Digitizer.h index d5e2a3fb9383f..0fdda0a6e8620 100644 --- a/Detectors/FIT/FV0/simulation/include/FV0Simulation/Digitizer.h +++ b/Detectors/FIT/FV0/simulation/include/FV0Simulation/Digitizer.h @@ -92,6 +92,7 @@ class Digitizer static Double_t PmtResponse(Double_t x); static Double_t PmtResponse(Double_t* x, Double_t*); static Double_t SinglePhESpectrum(Double_t* x, Double_t* par); + float getDistFromCellCenter(UInt_t cellId, double hitx, double hity); ClassDefNV(Digitizer, 1); }; diff --git a/Detectors/FIT/FV0/simulation/src/Detector.cxx b/Detectors/FIT/FV0/simulation/src/Detector.cxx index 09f6a535b0859..2a2ec90c637ec 100644 --- a/Detectors/FIT/FV0/simulation/src/Detector.cxx +++ b/Detectors/FIT/FV0/simulation/src/Detector.cxx @@ -267,7 +267,7 @@ void Detector::ConstructGeometry() { LOG(INFO) << "FV0: Constructing geometry"; createMaterials(); - mGeometry = new Geometry(Geometry::eFull); + mGeometry = Geometry::instance(Geometry::eFull); // mGeometry->enableComponent(Geometry::eScintillator, false); // mGeometry->enableComponent(Geometry::ePlastics, false); // mGeometry->enableComponent(Geometry::eFibers, false); diff --git a/Detectors/FIT/FV0/simulation/src/Digitizer.cxx b/Detectors/FIT/FV0/simulation/src/Digitizer.cxx index 1d2dd711eb3ba..0468384336570 100644 --- a/Detectors/FIT/FV0/simulation/src/Digitizer.cxx +++ b/Detectors/FIT/FV0/simulation/src/Digitizer.cxx @@ -9,6 +9,7 @@ // or submit itself to any jurisdiction. #include "FV0Simulation/Digitizer.h" +#include "FV0Base/Geometry.h" #include #include @@ -122,6 +123,9 @@ void Digitizer::process(const std::vector& hits) Float_t const t = hit.GetTime() * 1e9 + FV0DigParam::Instance().pmtTransitTime; Float_t const charge = TMath::Qe() * FV0DigParam::Instance().pmtGain * mBinSize / mPmtTimeIntegral; + // Example how to access the fields necessary to split the readout of cells in ring5 to two PMTs + // LOG(INFO) << detId << " " << Geometry::instance()->isRing5(detId) << " " << hit.GetX() << ", " << hit.GetY() << " " << getDistFromCellCenter(detId, hit.GetX(), hit.GetY()); + auto& analogSignal = mPmtChargeVsTime[detId]; for (Int_t iPhE = 0; iPhE < nPhE; ++iPhE) { @@ -254,3 +258,20 @@ Double_t Digitizer::SinglePhESpectrum(Double_t* x, Double_t*) return (TMath::Poisson(x[0], FV0DigParam::Instance().pmtNbOfSecElec) + FV0DigParam::Instance().pmtTransparency * TMath::Poisson(x[0], 1.0)); } + +// The Distance is positive for top half-sectors (when the hit position is above the cell center (has higher y)) +// TODO: performance check needed +float Digitizer::getDistFromCellCenter(UInt_t cellId, double hitx, double hity) +{ + Geometry* geo = Geometry::instance(); + + // Parametrize the line (ax+by+c=0) that crosses the detector center and the cell's middle point + Point3D* pCell = &geo->getCellCenter(cellId); + float x0, y0, z0; + geo->getGlobalPosition(x0, y0, z0); + double a = -(y0 - pCell->Y()) / (x0 - pCell->X()); + double b = 1; + double c = -(y0 - a * x0); + // Return the distance from hit to this line + return (a * hitx + b * hity + c) / TMath::Sqrt(a * a + b * b); +}