From ace9574fd6c40f0d85be7648a63dd7f0e7bedaa7 Mon Sep 17 00:00:00 2001
From: wiechula <Jens.Wiechula@cern.ch>
Date: Mon, 18 Nov 2024 23:03:06 +0100
Subject: [PATCH] Exclude electrons created between vessels and field strips

---
 .../base/include/TPCBase/ParameterDetector.h  |  9 +++----
 Detectors/TPC/simulation/src/Detector.cxx     | 24 +++++++++++++++++++
 2 files changed, 29 insertions(+), 4 deletions(-)

diff --git a/Detectors/TPC/base/include/TPCBase/ParameterDetector.h b/Detectors/TPC/base/include/TPCBase/ParameterDetector.h
index 7510d7246acd5..2762f6ff67d31 100644
--- a/Detectors/TPC/base/include/TPCBase/ParameterDetector.h
+++ b/Detectors/TPC/base/include/TPCBase/ParameterDetector.h
@@ -28,11 +28,12 @@ namespace tpc
 
 struct ParameterDetector : public o2::conf::ConfigurableParamHelper<ParameterDetector> {
 
-  float TPClength = 250.f;     ///< Length of the TPC [cm]
+  float TPClength = 250.f;       ///< Length of the TPC [cm]
   float TPCRecoWindowSim = 1.5f; ///< length of the reconstruction window in units of drift time of the TPC in simulation (Neutron capture process can extend up to 30-40 TPC drift time)
-  float PadCapacitance = 0.1f; ///< Capacitance of a single pad [pF]
-  TimeBin TmaxTriggered = 550; ///< Maximum time bin in case of triggered readout mode
-  float DriftTimeOffset = 7.3; ///< drift time offset in time bins (we observe ~2.4\mus before October 2023 and ~1.45 \mus after)
+  float PadCapacitance = 0.1f;   ///< Capacitance of a single pad [pF]
+  TimeBin TmaxTriggered = 550;   ///< Maximum time bin in case of triggered readout mode
+  float DriftTimeOffset = 7.3;   ///< drift time offset in time bins (we observe ~2.4\mus before October 2023 and ~1.45 \mus after)
+  bool ExcludeFCGap = true;      ///< exclude electrons created in the gap between the IFC vessel and OFC vessel and FC strips
 
   O2ParamDef(ParameterDetector, "TPCDetParam");
 };
diff --git a/Detectors/TPC/simulation/src/Detector.cxx b/Detectors/TPC/simulation/src/Detector.cxx
index 367abdc1a753a..e261424c41332 100644
--- a/Detectors/TPC/simulation/src/Detector.cxx
+++ b/Detectors/TPC/simulation/src/Detector.cxx
@@ -13,6 +13,7 @@
 #include "TPCSimulation/Detector.h"
 #include "TPCSimulation/Point.h"
 #include "TPCBase/ParameterGas.h"
+#include "TPCBase/ParameterDetector.h"
 
 #include "DetectorsBase/Stack.h"
 #include "SimulationDataFormat/TrackReference.h"
@@ -104,6 +105,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
 {
   mStepCounter++;
   auto& gasParam = ParameterGas::Instance();
+  auto& detParam = ParameterDetector::Instance();
   const Int_t kMaxDistRef = 15;       // maximal difference between 2 stored references - the parameter should be 15 cm as default
   static Double_t lastReferenceR = 0; // keeps last reference point in radius (cm)
 
@@ -140,6 +142,28 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
   // TODO: Temporary hack to process only one sector
   // if (sectorID != 0) return kFALSE;
 
+  // ---| remove clusters between the IFC and the FC strips |---
+  // those should not enter the active readout area
+  // do coarse selection before, to limit number of transformations
+  if (detParam.ExcludeFCGap) {
+    const auto rCluster = std::sqrt(position.X() * position.X() + position.Y() * position.Y());
+    const float rodRin = 81.5 + 2.2;    // radial position of the inner field cage rods + radial size of the field cage rods
+    const float rodRout = 254.25 + 2.2; // radial position of the outer field cage rods + radial size of the field cage rods
+    const float fcLxIn = 82.428409;     // position of the inner FC strips in local x = cos(10 deg) * rodRin;
+    const float fcLxOut = 252.55395;    // position of the outer FC strips in local x = cos(10 deg) * rodRin;
+
+    if (rCluster < rodRin || rCluster > fcLxOut) {
+      const int sectorIDnonShift = static_cast<int>(Sector::ToSector(position.X(), position.Y(), position.Z()));
+      const double alpha = TMath::DegToRad() * (10. + sectorIDnonShift * 20.);
+      const double cs = std::cos(-alpha), sn = std::sin(-alpha);
+      const auto localX = position.X() * cs - position.Y() * sn;
+      // fine cut
+      if (localX < fcLxIn || localX > fcLxOut) {
+        return kFALSE;
+      }
+    }
+  }
+
   // ---| momentum and beta gamma |---
   static TLorentzVector momentum; // static to make avoid creation/deletion of this expensive object
   fMC->TrackMomentum(momentum);