Add geometrical TPC protection against misfired loopers

Author: Marco Giacalone

Generators/include/Generators/TPCLoopers.h

@@ -107,8 +107,16 @@ class GenTPCLoopers
 
   void SetAdjust(float adjust = 0.f);
 
+  void setGeomProtection(bool protect);
+
+  // check if a vertex lies in the TPC volume where ionisation can be recorded
+  bool isInTPCActiveVolume(double vx, double vy) const;
+
   unsigned int getNLoopers() const { return (mNLoopersPairs + mNLoopersCompton); }
 
+  // loopers dropped by the geometrical protection since the last reset
+  unsigned int getNSkipped() const { return mNSkippedPairs + mNSkippedCompton; }
+
  private:
   std::unique_ptr<ONNXGenerator> mONNX_pair = nullptr;
   std::unique_ptr<ONNXGenerator> mONNX_compton = nullptr;
@@ -137,6 +145,9 @@ class GenTPCLoopers
   double mTimeEnd = 0.0;                                          // Time limit for the last event
   float mLoopsFractionPairs = 0.08;                               // Fraction of loopers from Pairs
   int mInteractionRate = 50000;                                   // Interaction rate in Hz
+  bool mGeomProtection = true;                                    // Skip loopers generated outside the TPC active volume
+  unsigned int mNSkippedPairs = 0;                                // Pairs dropped by the geometrical protection
+  unsigned int mNSkippedCompton = 0;                              // Compton electrons dropped by the geometrical protection
 };
 #endif // GENERATORS_WITH_TPCLOOPERS
 

Generators/include/Generators/TPCLoopersParam.h

@@ -45,6 +45,7 @@ struct GenTPCLoopersParam : public o2::conf::ConfigurableParamHelper<GenTPCLoope
   float multiplier[2] = {1., 1.};                                                                 // multiplier for pairs and compton loopers for Poissonian and Gaussian sampling
   unsigned int fixedNLoopers[2] = {1, 1};                                                         // fixed number of loopers coming from pairs and compton electrons - valid if flat gas is false and both Poisson and Gaussian params files are empty
   float adjust_flatgas = 0.f;                                                                     // adjustment for the number of flat gas loopers per orbit (in percentage, e.g. -0.1 = -10%) [-1, inf)]
+  bool geomProtection = true;                                                                     // skip loopers generated outside the TPC active volume
   O2ParamDef(GenTPCLoopersParam, "GenTPCLoopers");
 };
 

Generators/src/Generator.cxx

@@ -195,6 +195,7 @@ bool Generator::initTPCLoopersGen()
   try {
     // Create the TPC loopers generator with the provided parameters
     mTPCLoopersGen = new o2::eventgen::GenTPCLoopers(model_pairs, model_compton, poisson, gauss, scaler_pair, scaler_compton);
+    mTPCLoopersGen->setGeomProtection(loopersParam.geomProtection);
     const auto& intrate = loopersParam.intrate;
     // Configure the generator with flat gas loopers defined per orbit with clusters/track info
     // If intrate is negative (default), automatic IR from collisioncontext.root will be used
@@ -260,6 +261,10 @@ Bool_t
     mParticles.insert(mParticles.end(), looperParticles.begin(), looperParticles.end());
 
     LOG(debug) << "Added " << looperParticles.size() << " looper particles";
+    const auto skippedLoopers = mTPCLoopersGen->getNSkipped();
+    if (skippedLoopers > 0) {
+      LOG(debug) << "Geometrical protection skipped " << skippedLoopers << " loopers outside the TPC active volume";
+    }
   }
 #endif
   return kTRUE;

Generators/src/TPCLoopers.cxx

@@ -126,6 +126,40 @@ namespace o2
 namespace eventgen
 {
 
+namespace
+{
+// Radial limits of the region from which a looper can still reach the TPC
+// sensitive gas. The field cage positions are those used by the
+// "ExcludeFCGap" selection in o2::tpc::Detector::ProcessHits()
+// A looper can enter the sensitive gas from just outside it, so an additional margin is set
+// with a factor two over the largest radial excursion which was observed in validation (~4.2 cm).
+//
+// No cut is applied on z because loopers spiral the field lines and a vertex as far as |z| = 283 cm
+// feeds hits into the gas. A cut here would discard loopers that produce TPC signals.
+constexpr double kFcLxIn = 82.428409;         // cm, inner field cage strips
+constexpr double kRodROut = 254.25 + 2.2;     // cm, outer field cage rods plus their radial size
+constexpr double kLooperRadialReach = 10.;    // cm, margin for the helix sweep
+constexpr double kTPCActiveRMin = kFcLxIn - kLooperRadialReach;
+constexpr double kTPCActiveRMax = kRodROut + kLooperRadialReach;
+} // namespace
+
+bool GenTPCLoopers::isInTPCActiveVolume(double vx, double vy) const
+{
+  const double vt = std::sqrt(vx * vx + vy * vy);
+  return (vt >= kTPCActiveRMin && vt <= kTPCActiveRMax);
+}
+
+void GenTPCLoopers::setGeomProtection(bool protect)
+{
+  mGeomProtection = protect;
+  if (mGeomProtection) {
+    LOG(debug) << "TPC loopers geometrical protection: ON (accepting vertices with "
+              << kTPCActiveRMin << " <= Vt <= " << kTPCActiveRMax << " cm)";
+  } else {
+    LOG(warning) << "TPC loopers geometrical protection: OFF - loopers will be generated outside the TPC active volume as well.";
+  }
+}
+
 GenTPCLoopers::GenTPCLoopers(std::string model_pairs, std::string model_compton,
                              std::string poisson, std::string gauss, std::string scaler_pair,
                              std::string scaler_compton)
@@ -267,11 +301,20 @@ std::vector<TParticle> GenTPCLoopers::importParticles()
   std::vector<TParticle> particles;
   const double mass_e = TDatabasePDG::Instance()->GetParticle(11)->Mass();
   const double mass_p = TDatabasePDG::Instance()->GetParticle(-11)->Mass();
+  mNSkippedPairs = 0;
+  mNSkippedCompton = 0;
   // Get looper pairs from the event
   for (auto& pair : mGenPairs) {
     double px_e, py_e, pz_e, px_p, py_p, pz_p;
     double vx, vy, vz, time;
     double e_etot, p_etot;
+    // The generative model is not currently fully constrained to the TPC geometry, so it places
+    // significant fraction of the vertices outside the drift gas.
+    // These are now dropped before they reach the transport.
+    if (mGeomProtection && !isInTPCActiveVolume(pair[6], pair[7])) {
+      mNSkippedPairs++;
+      continue;
+    }
     px_e = pair[0];
     py_e = pair[1];
     pz_e = pair[2];
@@ -301,6 +344,10 @@ std::vector<TParticle> GenTPCLoopers::importParticles()
     double px, py, pz;
     double vx, vy, vz, time;
     double etot;
+    if (mGeomProtection && !isInTPCActiveVolume(compton[3], compton[4])) {
+      mNSkippedCompton++;
+      continue;
+    }
     px = compton[0];
     py = compton[1];
     pz = compton[2];