Merge 89f314b into sapling-pr-archive-ktf

Author: ktf

Detectors/Passive/src/FrameStructure.cxx

@@ -1789,7 +1789,21 @@ void FrameStructure::ConstructGeometry()
   ppgon[9] = ppgon[6];
 
   vmc->Gsvolu("BBMO", "PGON", kAir, ppgon, 10);
-  vmc->Gsdvn("BBCE", "BBMO", 18, 2);
+
+  // The 18 sectors, placed one by one rather than made with a phi division.
+  // Geant4 has no faithful representation of a phi division of a polyhedra: it
+  // divides by the number of sides and ignores the requested width and offset,
+  // so the sector contents end up half a sector away from where TGeo puts them.
+  const int kNSectors = 18;
+  const float kSectorDphi = 360. / kNSectors;
+  TGeoPgon* shBBCE = new TGeoPgon(-kSectorDphi / 2., kSectorDphi, 1, 2);
+  shBBCE->DefineSection(0, -kBBMdz / 2., kBBMRin, kBBMRou);
+  shBBCE->DefineSection(1, kBBMdz / 2., kBBMRin, kBBMRou);
+  TGeoVolume* voBBCE = new TGeoVolume("BBCE", shBBCE, kMedAir);
+  TGeoVolume* voBBMO = gGeoManager->GetVolume("BBMO");
+  for (i = 0; i < kNSectors; i++) {
+    voBBMO->AddNode(voBBCE, i + 1, new TGeoRotation("", (i + 0.5) * kSectorDphi, 0., 0.));
+  }
 
   // CBL ////////////////////////////////////////////////////////
   //

Detectors/Upgrades/ALICE3/IOTOF/macros/CheckDigitsIOTOF.C

@@ -75,7 +75,8 @@ void addTLines(float pitch)
   gPad->Update();
 }
 
-void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root")
+void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfile = "o2sim_HitsTF3.root", std::string inputGeom = "o2sim_geometry.root",
+                      std::string cfgStr = "IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false;")
 {
   gStyle->SetPalette(55);
 
@@ -85,7 +86,7 @@ void CheckDigitsIOTOF(std::string digifile = "tf3digits.root", std::string hitfi
   using o2::iotof::Digit;
   using o2::itsmft::Hit;
 
-  o2::conf::ConfigurableParam::updateFromString("IOTOFBase.segmentedInnerTOF=true;IOTOFBase.segmentedOuterTOF=true;IOTOFBase.enableForwardTOF=false;IOTOFBase.enableBackwardTOF=false");
+  o2::conf::ConfigurableParam::updateFromString(cfgStr);
 
   auto seg = o2::iotof::Segmentation::Instance();
 

Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/DPLDigitizerParam.h

@@ -34,7 +34,7 @@ struct DPLDigitizerParam : public o2::conf::ConfigurableParamHelper<DPLDigitizer
                                           ///< The efficiency map is currently available at /alice/cern.ch/user/g/glucia/ALICE3/IOTOF/pixelEfficiency/PixelEfficiencyMap_TH2.root. FIXME to be removed once switch to CCDBFetcher
   int chargeThreshold = 100;              ///< charge threshold in Nelectrons
   int minChargeToAccount = 7;             ///< minimum charge contribution to account
-  int nSimSteps = 1;                      ///< number of steps in response simulation
+  int nSimSteps = 10;                     ///< number of steps in response simulation
   float energyToNElectrons = 1. / 3.6e-9; // conversion of eloss to Nelectrons
   int responseMatrixSize = 1;             ///< size of the response matrix (odd number)
 

Detectors/Upgrades/ALICE3/IOTOF/simulation/include/IOTOFSimulation/Digitizer.h

@@ -83,7 +83,7 @@ class Digitizer : public TObject
   void registerDigits(Chip& chip, uint32_t roFrame, double time, int nROF,
                       uint16_t row, uint16_t col, int nElectrons, o2::MCCompLabel& label);
 
-  void stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& rowStart, int& colStart, int& rowSpan, int& colSpan);
+  void stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float**& avgHitLocalX, float**& avgHitLocalZ, int& rowStart, int& colStart, int& rowSpan, int& colSpan);
 
   /// Apply time smearing to simulate detector resolution
   double smearTime(double time) const;

Detectors/Upgrades/ALICE3/IOTOF/simulation/src/Digitizer.cxx

@@ -104,32 +104,18 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
 
   // Get detector element ID
   const int chipID = hit.GetDetectorID();
+  if (chipID < 0 || chipID >= mGeometry->getSize() || mGeometry->getSize() < 1) {
+    LOG(debug) << "Invalid detector ID: " << chipID << ", geometry size: " << mGeometry->getSize();
+    return; // invalid detector ID
+  }
+  const int subdetectorID = mGeometry->getIOTOFLayer(chipID);
+
   auto& chip = mChips[chipID];
   if (chip.isDisabled()) {
     LOG(debug) << "Hit rejected because chip " << chipID << " is disabled";
     return;
   }
 
-  // middle position of the hit in the sensor frame
-  const auto& matrix = mGeometry->getMatrixL2G(chipID);
-  auto xyzPositionStart = matrix ^ hit.GetPosStart();
-  auto xyzPositionEnd = matrix ^ hit.GetPos();
-  const auto xMid = 0.5f * (xyzPositionStart.X() + xyzPositionEnd.X());
-  const auto zMid = 0.5f * (xyzPositionStart.Z() + xyzPositionEnd.Z());
-  // move this to the local pixel coordinates for the efficiency map
-  int row, col;
-  float xPixelCenter, zPixelCenter;
-  if (!sSegmentation->localToDetector(xMid, zMid, row, col, mGeometry->getIOTOFLayer(chipID))) {
-    LOG(debug) << "Hit rejected because position (" << xMid << ", " << zMid << ") is outside the active area of chip " << chipID;
-    return; // hit is outside the active area
-  }
-  sSegmentation->detectorToLocalUnchecked(row, col, xPixelCenter, zPixelCenter, mGeometry->getIOTOFLayer(chipID));
-
-  if (!isEfficient(xMid - xPixelCenter, zMid - zPixelCenter)) {
-    LOG(debug) << "Hit rejected by efficiency cut";
-    return;
-  }
-
   // Convert energy loss to charge (number of electrons)
   float energyLoss = hit.GetEnergyLoss(); // in GeV
   int charge = energyToCharge(energyLoss);
@@ -149,20 +135,18 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
   double hitTimeWrtBC = hitTime + eventTimeInBC;          // hit time wrt bc
   double smearedTime = smearTime(hitTimeWrtBC);
 
-  if (chipID < 0 || chipID >= mGeometry->getSize() || mGeometry->getSize() < 1) {
-    LOG(debug) << "Invalid detector ID: " << chipID << ", geometry size: " << mGeometry->getSize();
-    return; // invalid detector ID
-  }
-
   // Create the digit with time information
   o2::MCCompLabel label(hit.GetTrackID(), evID, srcID, false);
   const int roFrameAbs = 0; // For now, we can set this to 0 or calculate based on time if needed
   const int nROF = 1;       // For now, we can assume the signal is contained in one ROF, this can be extended to multiple ROFs based on the time
 
   float** respMatrix = nullptr;
+  float** avgHitLocalX = nullptr;
+  float** avgHitLocalZ = nullptr;
   int rowStart = 0, colStart = 0, rowSpan = 0, colSpan = 0;
-  stepping(hit, respMatrix, rowStart, colStart, rowSpan, colSpan);
+  stepping(hit, respMatrix, avgHitLocalX, avgHitLocalZ, rowStart, colStart, rowSpan, colSpan);
 
+  float xPixelCenter = 0.0f, zPixelCenter = 0.0f;
   for (int irow = rowSpan; irow--;) {
     uint16_t rowIS = irow + rowStart;
     for (int icol = colSpan; icol--;) {
@@ -171,6 +155,14 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
       if (!nEleResp) {
         continue;
       }
+
+      // Apply efficiency cut based on the hit segment mean position relative to the pixel center
+      sSegmentation->detectorToLocal(rowIS, colIS, xPixelCenter, zPixelCenter, subdetectorID);
+      if (!isEfficient(avgHitLocalX[irow][icol] - xPixelCenter, avgHitLocalZ[irow][icol] - zPixelCenter)) {
+        LOG(debug) << "Hit rejected by efficiency cut at pixel (" << rowIS << ", " << colIS << ") in chip " << chipID;
+        continue;
+      }
+
       const int nElectronsSampled = gRandom->Poisson(electronsPerStep * nEleResp);
       // Noise can be added here if needed
 
@@ -181,12 +173,17 @@ void Digitizer::processHit(const o2::itsmft::Hit& hit, int evID, int srcID)
 
   for (int irow = 0; irow < rowSpan; ++irow) {
     delete[] respMatrix[irow];
+    delete[] avgHitLocalX[irow];
+    delete[] avgHitLocalZ[irow];
   }
   delete[] respMatrix;
+  delete[] avgHitLocalX;
+  delete[] avgHitLocalZ;
 }
 
-void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& rowStart, int& colStart, int& rowSpan, int& colSpan)
+void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, float**& avgHitLocalX, float**& avgHitLocalZ, int& rowStart, int& colStart, int& rowSpan, int& colSpan)
 {
+  LOG(debug) << "\n\nPerforming stepping";
   const int chipID = hit.GetDetectorID();
   const auto& matrix = mGeometry->getMatrixL2G(chipID);
   const int subdetectorID = mGeometry->getIOTOFLayer(chipID);
@@ -238,27 +235,52 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& r
   colSpan = colEnd - colStart + 1;
 
   respMatrix = new float*[rowSpan];
+  avgHitLocalX = new float*[rowSpan];
+  avgHitLocalZ = new float*[rowSpan];
   for (int i = 0; i < rowSpan; ++i) {
     respMatrix[i] = new float[colSpan]();
+    avgHitLocalX[i] = new float[colSpan]();
+    avgHitLocalZ[i] = new float[colSpan]();
   }
 
-  int rowPrev = -1, colPrev = -1, row = 0, col = 0;
-  if (!respMatrix || rowSpan <= 0 || colSpan <= 0) {
+  if (!respMatrix || !avgHitLocalX || !avgHitLocalZ || rowSpan <= 0 || colSpan <= 0) {
     return;
   }
   if (nSkip) {
     nSteps -= nSkip;
   }
 
-  auto& currentPosLocal = xyzPositionStart;
+  int rowPrev = -1, colPrev = -1, row = 0, col = 0;
+  auto pixelCurrentPosLocal = xyzPositionStart;
+  auto pixelStartPosLocal = xyzPositionStart;
   for (int iStep = nSteps; iStep--;) {
-    sSegmentation->localToDetector(currentPosLocal.X(), currentPosLocal.Z(), row, col, subdetectorID);
+
+    // Step does not contribute if it is in the passive area
+    if (!sSegmentation->localToDetector(pixelCurrentPosLocal.X(), pixelCurrentPosLocal.Z(), row, col, subdetectorID)) {
+      LOG(debug) << "Step is in passive area: (" << pixelCurrentPosLocal.X() << ", " << pixelCurrentPosLocal.Z() << ") is outside the active area of chip " << subdetectorID;
+      pixelCurrentPosLocal += stepVector;
+      continue;
+    }
+
+    // The step has reached another pixel, compute mean hit segment positions
+    // for pixel efficiency evaluation and reset the start position for the next pixel
     if (row != rowPrev || col != colPrev) {
+
+      // Finalize the previous pixel
+      if (rowPrev != -1 && colPrev != -1) {
+        const int irow = rowPrev - rowStart;
+        const int icol = colPrev - colStart;
+        avgHitLocalX[irow][icol] = 0.5f * (pixelStartPosLocal.X() + pixelCurrentPosLocal.X() - stepVector.X());
+        avgHitLocalZ[irow][icol] = 0.5f * (pixelStartPosLocal.Z() + pixelCurrentPosLocal.Z() - stepVector.Z());
+      }
+
+      // Start the new pixel
       rowPrev = row;
       colPrev = col;
+      pixelStartPosLocal = pixelCurrentPosLocal;
     }
 
-    currentPosLocal += stepVector; // Move to the next step position
+    pixelCurrentPosLocal += stepVector; // Move to the next step position
 
     for (int irow = digitizerParams.responseMatrixSize; irow--;) {
       int rowDest = row + irow - (digitizerParams.responseMatrixSize / 2) - rowStart; // destination row in the respMatrix
@@ -274,6 +296,14 @@ void Digitizer::stepping(const o2::itsmft::Hit& hit, float**& respMatrix, int& r
       }
     }
   }
+
+  // Finalize the last pixel
+  if (rowPrev != -1 && colPrev != -1) {
+    const int irow = rowPrev - rowStart;
+    const int icol = colPrev - colStart;
+    avgHitLocalX[irow][icol] = 0.5f * (pixelStartPosLocal.X() + pixelCurrentPosLocal.X() - stepVector.X());
+    avgHitLocalZ[irow][icol] = 0.5f * (pixelStartPosLocal.Z() + pixelCurrentPosLocal.Z() - stepVector.Z());
+  }
 }
 
 //_______________________________________________________________________
@@ -407,7 +437,7 @@ void Digitizer::registerDigits(Chip& chip, uint32_t roFrame, double time, int nR
   int tdc = int((time - nbc * o2::constants::lhc::LHCBunchSpacingNS) / digitizerParams.tdcBin);
   nbc += mEventTime.toLong();
 
-  LOG(debug) << nbc << "\t" << tdc;
+  LOG(debug) << "nbc: " << nbc << "\ttdc: " << tdc;
   double absoluteTime = tdc * digitizerParams.tdcBin * 1.e-9 + nbc * o2::constants::lhc::LHCBunchSpacingNS;
 
   auto key = o2::iotof::Digit::getOrderingKey(nbc, row, col);

Framework/AnalysisSupport/src/AODJAlienReaderHelpers.cxx

@@ -32,6 +32,7 @@
 #include "Framework/EndOfStreamContext.h"
 #include "Framework/DeviceSpec.h"
 #include "Framework/RawDeviceService.h"
+#include "Framework/RuntimeError.h"
 #include "Framework/DataSpecUtils.h"
 #include "Framework/MessageContext.h"
 #include "Framework/Signpost.h"
@@ -131,7 +132,9 @@ static std::string describeException(std::exception const& exception)
   try {
     std::rethrow_if_nested(exception);
   } catch (std::exception const& nested) {
-    description += ": " + describeException(nested);
+    description += fmt::format(": {}", describeException(nested));
+  } catch (RuntimeErrorRef const& ref) {
+    description += fmt::format(": {}", error_from_ref(ref).what);
   } catch (...) {
     description += ": unknown exception";
   }
@@ -275,6 +278,7 @@ AlgorithmSpec AODJAlienReaderHelpers::rootFileReaderCallback(ConfigContext const
         auto skippedTimeframes = ++totalInvalidReadSkipped;
         LOGP(error, "Invalid AOD read for table {}: fileCounter {}, timeFrame {}. Skipping timeframe (skipped timeframes: {}). Reason: {}",
              concrete.origin.as<std::string>(), fcnt, ntf, skippedTimeframes, describeException(e));
+        clean_all_runtime_errors();
         didir->markTimeFrameSkipped(header::DataHeader(concrete.description, concrete.origin, concrete.subSpec), ntf);
         arrowContext.clear();
         messageContext.discard();

Framework/AnalysisSupport/src/DataInputDirector.cxx

@@ -574,6 +574,9 @@ bool DataInputDescriptor::readTree(DataAllocator& outputs, header::DataHeader dh
   }
 
   auto schemaOpt = format->Inspect(fullpath);
+  if (!schemaOpt.ok()) {
+    throw InvalidAODReadError(fmt::format("Unable to inspect tree {}: {}", treename, schemaOpt.status().ToString()));
+  }
   auto physicalSchema = schemaOpt;
   std::vector<std::shared_ptr<arrow::Field>> fields;
   for (auto& original : (*schemaOpt)->fields()) {