[ITS3] Update default parameters and adapt code to work with fourth layer (#11227)

* Change default parameters

* Avoid to get GRPECS from CCDB

* Fix simulation and digitisation with 4th later

* Adapt clusterer to work with 4th layer

* Adapt ITS3 tracker to work with 4th layer

* Update README

* Fix typo

* Please consider the following formatting changes

* Move modifications to ITS code in dedicated PR

* Remove leftover

* Fix conflicts

* Fix conflicts

---------

Co-authored-by: ALICE Action Bot <alibuild@cern.ch>

Author: fgrosa

Detectors/Upgrades/ITS3/README.md

@@ -13,7 +13,7 @@ Events can be simulated using the `o2-sim` workflow. To include ITS3 in the simu
 
 The following command can be used to generate heavy-ion collisions:
 ```bash
-o2-sim -j 1  \
+o2-sim -j 8  \
 -n 10 -g pythia8hi --field ccdb \
 --configKeyValues "Diamond.width[2]=6.;DescriptorInnerBarrelITS3.mVersion=ThreeLayers" \
 --run 311935
@@ -22,14 +22,16 @@ In the previous command:
 - `-j` is used to set the number of threads;
 - `-n` is used to set the number of events to simulate;
 - `-g` is used to set the event generator, in this case `pythia8hi`. To simulate pp collisions one can use `pythia8pp`.
-- `--configKeyValues` is needed to set internal parameters of the workflow. Among these parameters, the geometry of the ITS3 inner barrel can be set.
+- `--configKeyValues` is needed to set internal parameters of the workflow. Among these parameters, the geometry of the ITS3 inner barrel can be set:
+  - `DescriptorInnerBarrel.mVersion` is the geometry version of the ITS3 inner barrel (`ThreeLayersNoDeadZones`, `ThreeLayers`, or `FourLayers`)
+  - `DescriptorInnerBarrel.mRadii` is a 4-element vector with the radii of the ITS3 layers
+  - `DescriptorInnerBarrel.mLength` is the length of the ITS3 in the Z direction
+  - `DescriptorInnerBarrel.mGapY` is a 4-element vector with the values of gap between the two hemicylinders, described by a translation of the hemicylinders in the vertical direction
+  - `DescriptorInnerBarrel.mGapPhi` is a 4-element vector with the values of gap of azimuthal angle between the two hemicylinders, described by the maximum distance between the two hemicylinders. Differently from `mGapY`, in this case there no shift in the vertical direction, but a smaller coverage in the azimuthal angle of the half layers.
+  - `DescriptorInnerBarrel.mGapXDirection4thLayer` is the gap in the horizontal direction for the fourth layer, analogous to the `mGapY`.
+  - `SuperAlpideParams.mDetectorThickness` is the thickness of the chip
 - `--run` is needed to set the run number.
 
-Currently, three different geometries of the ITS3 inner barrel are available:
-- `ThreeLayersNoDeadZones`
-- `ThreeLayers`
-- `FourLayers`
-
 The run number is needed to retrieve objects from the CCDB. There are specific ranges of run-numbers, according to the collision system and to the selected geometry if the ITS3 inner barrel:
 
 - **pp** collisions:
@@ -42,7 +44,6 @@ The run number is needed to retrieve objects from the CCDB. There are specific r
   - 311934—311966 (`ThreeLayers`)
   - 311967—311999 (`FourLayers`)
 
-
 ### Using external generators based on AliRoot
 It is also possible to simulate heavy-ion collision using external generators based on AliRoot. In this case, it is necessary to load both O2 and AliROOT (the order is important):
 
@@ -77,17 +78,28 @@ o2-sim-digitizer-workflow -b --run --interactionRate 50000 \
 ```
 As above, it is important to set the correct interaction rate and run number.
 
+In addition, some parameters related to the segmentation of the chips can be set with the `--configKeyValues` argument:
+- `SuperAlpideParams.mPitchCol` is the pitch of the column (Z direction)
+- `SuperAlpideParams.mPitchRow` is the pitch of the row (r$\phi$ direction)
+- `SuperAlpideParams.mDetectorThickness` is the thickness of the chip
+
+If some parameters of the geometry changed in the simulation, they should be set similarly for the digitisation. In particular, the `mVersion`, `mRadii`, and `mLength` must be set if different from the default ones.
+
 # Reconstruction
 
 In this step, clustering, vertexing and tracking are performed. This is obtained with the `o2-its3-reco-workflow`:
 
 ```bash
 o2-its3-reco-workflow --tracking-mode async -b --run \
 -—configKeyValues "
-HBFUtils.runNumber=311935;ITSCATrackerParam.trackletsPerClusterLimit=20;ITSCATrackerParam.cellsPerC
-lusterLimit=20;ITSVertexerParam.lowMultXYcut2=0."
+HBFUtils.runNumber=311935;ITSCATrackerParam.trackletsPerClusterLimit=20;ITSCATrackerParam.cellsPerClusterLimit=20;ITSVertexerParam.lowMultXYcut2=0."
 ```
 
-As above, it is important to provide the correct run number using `-—configKeyValues`, to retrieve the correct files from the CCDB. The other internal parameters provided `-—configKeyValues` via are specific to the cased here considered (Pb-Pb) and are inherited from ITS2.
+As above, it is important to provide the correct run number using `-—configKeyValues`, to retrieve the correct files from the CCDB. The other internal parameters for the vertexer and the tracker are provided `-—configKeyValues` via are specific to the cased here considered (Pb-Pb) and are inherited from ITS2.
+
+If the `FourLayers` geometry was used in the simulation, it should be set also for the reconstruction to set properly the tracker to work with the additional layer. If something else of the geometry was set differently (`mRadii`, `mLength`, `mGapY`, `mGapPhi`, `mGapXDirection4thLayer`, or `mDetectorThickness`), it is necessary to remap the file with the geometry to replace the one on the CCDB, which would be different. This can be done by copying the `o2sim_geometry-aligned.root` file created during the simulation to a directory called `GLO/Config/GeometryAligned`, with the name `snapshot.root` in a local path of choice. Then, the following argument has to be added to the reco workflow: `--condition-remap "file://local_path=GLO/Config/GeometryAligned"`.
 
-> **_NOTE:_**  reconstruction for the `FourLayers` geometry is not implemented yet.
+> **_NOTE:_**  in order to make the reconstruction for the `FourLayers` geometry work, the following line has to be added to the CMake options in the [o2.sh](https://github.com/alisw/alidist/blob/master/o2.sh) recipe in alidist:
+```bash
+${ENABLE_ITS3_4L:+-DENABLE_ITS3_4L=ON}
+```

Detectors/Upgrades/ITS3/base/include/ITS3Base/DescriptorInnerBarrelITS3Param.h

@@ -39,11 +39,11 @@ enum class ITS3Version {
 struct DescriptorInnerBarrelITS3Param : public o2::conf::ConfigurableParamHelper<DescriptorInnerBarrelITS3Param> {
   ITS3Version mVersion = ITS3Version::None;
   int mBuildLevel{0};
-  double mGapY[4] = {0.1f, 0.1f, 0.1f, 0.1f};
-  double mGapPhi[4] = {0.f, 0.f, 0.f, 0.f};
+  double mGapY[4] = {0.f, 0.f, 0.f, 0.f};
+  double mGapPhi[4] = {0.1f, 0.1f, 0.1f, 0.1f};
   double mRadii[4] = {1.8f, 2.4f, 3.0f, 6.0f};
-  double mLength{27.f};
-  double mGapXDirection4thLayer{0.05f};
+  double mLength{26.f};
+  double mGapXDirection4thLayer{0.f};
   std::string const& getITS3LayerConfigString() const;
   O2ParamDef(DescriptorInnerBarrelITS3Param, "DescriptorInnerBarrelITS3");
 };

Detectors/Upgrades/ITS3/base/include/ITS3Base/SegmentationSuperAlpide.h

@@ -40,20 +40,20 @@ class SegmentationSuperAlpide
   }
   SegmentationSuperAlpide(int layer = 0) : SegmentationSuperAlpide(layer, SuperAlpideParams::Instance().mPitchCol, SuperAlpideParams::Instance().mPitchRow, SuperAlpideParams::Instance().mDetectorThickness, DescriptorInnerBarrelITS3Param::Instance().mLength, DescriptorInnerBarrelITS3Param::Instance().mRadii) {}
 
-  double mRadii[4] = {1.8f, 2.4f, 3.0f, 6.0f};                                                                                                                                    ///< radii for different layers
-  const double mLength;                                                                                                                                                           ///< chip length
-  const int mLayer;                                                                                                                                                               ///< chip layer
-  const float mPitchCol;                                                                                                                                                          ///< pixel column size
-  const float mPitchRow;                                                                                                                                                          ///< pixel row size
-  const float mDetectorLayerThickness;                                                                                                                                            ///< detector thickness
-  const int mNCols{static_cast<int>(std::ceil(mLength / mPitchCol))};                                                                                                             ///< number of columns
-  const int mNRows{static_cast<int>(std::ceil(double(mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThickness / 2.) * double(constants::math::PI) / double(mPitchRow)))}; ///< number of rows
-  const int mNPixels{mNRows * mNCols};                                                                                                                                            ///< total number of pixels
-  static constexpr float mPassiveEdgeReadOut = 0.;                                                                                                                                ///< width of the readout edge (Passive bottom)
-  static constexpr float mPassiveEdgeTop = 0.;                                                                                                                                    ///< Passive area on top
-  static constexpr float mPassiveEdgeSide = 0.;                                                                                                                                   ///< width of Passive area on left/right of the sensor
-  const float mActiveMatrixSizeCols{mPitchCol * mNCols};                                                                                                                          ///< Active size along columns
-  const float mActiveMatrixSizeRows{mPitchRow * mNRows};                                                                                                                          ///< Active size along rows
+  double mRadii[4] = {1.8f, 2.4f, 3.0f, 6.0f};                                                                                                                                                               ///< radii for different layers
+  const double mLength;                                                                                                                                                                                      ///< chip length
+  const int mLayer;                                                                                                                                                                                          ///< chip layer
+  const float mPitchCol;                                                                                                                                                                                     ///< pixel column size
+  const float mPitchRow;                                                                                                                                                                                     ///< pixel row size
+  const float mDetectorLayerThickness;                                                                                                                                                                       ///< detector thickness
+  const int mNCols{static_cast<int>(std::ceil(mLength / mPitchCol))};                                                                                                                                        ///< number of columns
+  const int mNRows{static_cast<int>(std::ceil(double(mRadii[mLayer] + mDetectorLayerThickness - mSensorLayerThickness / 2.) * double(constants::math::PI) / double(mPitchRow) * (mLayer == 3 ? 0.5 : 1.)))}; ///< number of rows
+  const int mNPixels{mNRows * mNCols};                                                                                                                                                                       ///< total number of pixels
+  static constexpr float mPassiveEdgeReadOut = 0.;                                                                                                                                                           ///< width of the readout edge (Passive bottom)
+  static constexpr float mPassiveEdgeTop = 0.;                                                                                                                                                               ///< Passive area on top
+  static constexpr float mPassiveEdgeSide = 0.;                                                                                                                                                              ///< width of Passive area on left/right of the sensor
+  const float mActiveMatrixSizeCols{mPitchCol * mNCols};                                                                                                                                                     ///< Active size along columns
+  const float mActiveMatrixSizeRows{mPitchRow * mNRows};                                                                                                                                                     ///< Active size along rows
 
   // effective thickness of sensitive layer, accounting for charge collection non-uniformity, https://alice.its.cern.ch/jira/browse/AOC-46
   static constexpr float mSensorLayerThicknessEff = 28.e-4;                                   ///< effective thickness of sensitive part

Detectors/Upgrades/ITS3/base/include/ITS3Base/SuperAlpideParams.h

@@ -24,8 +24,8 @@ namespace its3
 
 /// Segmentation parameters for Super ALPIDE chips
 struct SuperAlpideParams : public o2::conf::ConfigurableParamHelper<SuperAlpideParams> {
-  float mPitchCol = 20.e-4;          ///< Pixel column size (cm)
-  float mPitchRow = 20.e-4;          ///< Pixel row size (cm)
+  float mPitchCol = 29.24e-4;        ///< Pixel column size (cm) //FIXME: proxy value to get same resolution as ITS2 given incorrect sensor response
+  float mPitchRow = 26.88e-4;        ///< Pixel row size (cm) //FIXME: proxy value to get same resolution as ITS2 given incorrect sensor response
   float mDetectorThickness = 50.e-4; ///< Detector thickness (cm)
 
   // boilerplate

Detectors/Upgrades/ITS3/macros/test/CheckClustersITS3.C

@@ -38,8 +38,8 @@
 #include "DetectorsCommonDataFormats/DetectorNameConf.h"
 #endif
 
-void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.root", std::string hitfile = "o2sim_HitsIT3.root",
-                       std::string inputGeom = "", std::string dictfile = "", bool batch = true)
+void CheckClustersITS3(std::string clusfile = "o2clus_it3.root", std::string hitfile = "o2sim_HitsIT3.root",
+                       std::string inputGeom = "o2sim_geometry.root", std::string dictfile = "", bool batch = true)
 {
   gROOT->SetBatch(batch);
 
@@ -58,15 +58,6 @@ void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.r
   std::vector<HitVec*> hitVecPool;
   std::vector<MC2HITS_map> mc2hitVec;
 
-  // we assume that we have 2 chips per layer
-  const int nChipsPerLayer = 2;
-  std::vector<SegmentationSuperAlpide> segs{};
-  for (int iLayer{0}; iLayer < nITS3layers; ++iLayer) {
-    for (int iChip{0}; iChip < nChipsPerLayer; ++iChip) {
-      segs.push_back(SegmentationSuperAlpide(iLayer));
-    }
-  }
-
   const int QEDSourceID = 99; // Clusters from this MC source correspond to QED electrons
 
   TFile fout("CheckClusters.root", "recreate");
@@ -78,6 +69,13 @@ void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.r
   gman->fillMatrixCache(o2::math_utils::bit2Mask(o2::math_utils::TransformType::T2L, o2::math_utils::TransformType::T2GRot,
                                                  o2::math_utils::TransformType::L2G)); // request cached transforms
 
+  std::vector<SegmentationSuperAlpide> segs{};
+  for (int iLayer{0}; iLayer < gman->getNumberOfLayers() - 4; ++iLayer) {
+    for (int iChip{0}; iChip < gman->getNumberOfChipsPerLayer(iLayer); ++iChip) {
+      segs.push_back(SegmentationSuperAlpide(iLayer));
+    }
+  }
+
   // Hits
   TFile fileH(hitfile.data());
   TTree* hitTree = (TTree*)fileH.Get("o2sim");
@@ -183,13 +181,13 @@ void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.r
       auto chipID = cluster.getSensorID();
       if (pattID == o2::its3::CompCluster::InvalidPatternID || dict.isGroup(pattID)) {
         o2::itsmft::ClusterPattern patt(pattIt);
-        locC = dict.getClusterCoordinates(cluster, patt, false);
+        locC = dict.getClusterCoordinates(cluster, patt, false, segs.size());
         LOGP(info, "I am invalid and I am on chip {}", chipID);
       } else {
-        locC = dict.getClusterCoordinates(cluster);
+        locC = dict.getClusterCoordinates(cluster, segs.size());
         errX = dict.getErrX(pattID);
         errZ = dict.getErrZ(pattID);
-        if (chipID / nChipsPerLayer >= nITS3layers) {
+        if (chipID >= segs.size()) {
           errX *= Segmentation::PitchRow;
           errZ *= Segmentation::PitchCol;
         } else {
@@ -231,7 +229,7 @@ void CheckClustersITS3(int nITS3layers = 3, std::string clusfile = "o2clus_it3.r
       auto y0 = locHsta.Y(), dlty = locH.Y() - y0;
       auto z0 = locHsta.Z(), dltz = locH.Z() - z0;
 
-      if (chipID / nChipsPerLayer >= nITS3layers) {
+      if (chipID >= segs.size()) {
         auto r = (0.5 * (Segmentation::SensorLayerThickness - Segmentation::SensorLayerThicknessEff) - y0) / dlty;
         locH.SetXYZ(x0 + r * dltx, y0 + r * dlty, z0 + r * dltz);
       } else {

Detectors/Upgrades/ITS3/macros/test/CheckDigitsITS3.C

@@ -38,7 +38,7 @@
 
 #endif
 
-void CheckDigitsITS3(int nITS3layers = 3, std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "", std::string paramfile = "o2sim_par.root", bool batch = true)
+void CheckDigitsITS3(std::string digifile = "it3digits.root", std::string hitfile = "o2sim_HitsIT3.root", std::string inputGeom = "o2sim_geometry.root", std::string paramfile = "o2sim_par.root", bool batch = true)
 {
   gROOT->SetBatch(batch);
 
@@ -62,8 +62,8 @@ void CheckDigitsITS3(int nITS3layers = 3, std::string digifile = "it3digits.root
   const int nChipsPerLayer = 2;
 
   std::vector<SegmentationSuperAlpide> segs{};
-  for (int iLayer{0}; iLayer < nITS3layers; ++iLayer) {
-    for (int iChip{0}; iChip < nChipsPerLayer; ++iChip) {
+  for (int iLayer{0}; iLayer < gman->getNumberOfLayers() - 4; ++iLayer) {
+    for (int iChip{0}; iChip < gman->getNumberOfChipsPerLayer(iLayer); ++iChip) {
       segs.push_back(SegmentationSuperAlpide(iLayer));
     }
   }
@@ -172,7 +172,7 @@ void CheckDigitsITS3(int nITS3layers = 3, std::string digifile = "it3digits.root
 
       int chipID = (*digArr)[iDigit].getChipIndex();
 
-      if (chipID / nChipsPerLayer < nITS3layers) {
+      if (chipID < segs.size()) {
         float xFlat{0.f};
         segs[chipID].detectorToLocal(ix, iz, xFlat, z);
         segs[chipID].flatToCurved(xFlat, 0., x, y);
@@ -214,7 +214,7 @@ void CheckDigitsITS3(int nITS3layers = 3, std::string digifile = "it3digits.root
         int row, col;
         float xlc = 0., zlc = 0.;
 
-        if (chipID / nChipsPerLayer < nITS3layers) {
+        if (chipID < segs.size()) {
           segs[chipID].localToDetector(locH.X(), locH.Z(), row, col);
           segs[chipID].detectorToLocal(row, col, xlc, zlc);
         } else {

Detectors/Upgrades/ITS3/macros/test/CheckTracksITS3.C

@@ -73,6 +73,7 @@ void CheckTracksITS3(std::string tracfile = "o2trac_its3.root",
                      std::string clusfile = "o2clus_it3.root",
                      std::string kinefile = "o2sim_Kine.root",
                      std::string magfile = "o2sim_grp.root",
+                     std::string inputGeom = "o2sim_geometry.root",
                      bool batch = true)
 {
 
@@ -94,7 +95,7 @@ void CheckTracksITS3(std::string tracfile = "o2trac_its3.root",
   double orig[3] = {0., 0., 0.};
   float bz = field->getBz(orig);
   // Geometry
-  o2::base::GeometryManager::loadGeometry();
+  o2::base::GeometryManager::loadGeometry(inputGeom);
   auto gman = o2::its::GeometryTGeo::Instance();
 
   // MC tracks