AliceO2Group · ktf · Apr 28, 2020 · Apr 27, 2020 · Apr 27, 2020
@@ -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)
 
@@ -18,10 +18,12 @@
 #define ALICEO2_FV0_GEOMETRY_H_
 
 #include <vector>
+#include <array>
 
 #include <TGeoMatrix.h>
 #include <TGeoVolume.h>
 #include <TVirtualMC.h>
+#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<float>& getCellCenter(UInt_t cellId);
+  Point3D<float>& 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<std::string> 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<Point3D<float>, sNumberOfCells> mCellCenter;              ///< Center of each scintillator cell
+  std::array<Point3D<float>, sNumberOfReadoutChannels> mReadoutCenter; ///< Similar to mCellCenter, cells in r5 are additionally divided
+
+  static Geometry* sInstance; ///< \brief  Singleton instance
+
   ClassDefNV(Geometry, 1);
 };
 } // namespace fv0

@@ -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<float>& Geometry::getCellCenter(UInt_t cellId)
+{
+  return mCellCenter.at(cellId);
+}
+
+Point3D<float>& 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<float>* p = &mCellCenter.at(cellId);
+    p->SetCoordinates(x, y, sZGlobal);
+  }
+}
+
+void Geometry::initializeReadoutCenters()
+{
+  for (int channelId = 0; channelId < sNumberOfReadoutChannels; channelId++) {
+    Point3D<float>* 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;
+}
@@ -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);
 };

@@ -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);

@@ -9,6 +9,7 @@
 // or submit itself to any jurisdiction.
 
 #include "FV0Simulation/Digitizer.h"
+#include "FV0Base/Geometry.h"
 
 #include <TRandom.h>
 #include <algorithm>
@@ -122,6 +123,9 @@ void Digitizer::process(const std::vector<o2::fv0::Hit>& 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<float>* 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);
+}