Set automatic interaction rate from collision context

Author: jackal1-66

Generators/include/Generators/TPCLoopersParam.h

@@ -37,9 +37,9 @@ struct GenTPCLoopersParam : public o2::conf::ConfigurableParamHelper<GenTPCLoope
   std::string scaler_compton = "${O2_ROOT}/share/Generators/egconfig/ScalerComptonParams.json"; // file with scaler parameters for Compton scattering
   std::string nclxrate = "ccdb://Users/m/mgiacalo/ClustersTrackRatio";                          // file with clusters/rate information per orbit
   std::string colsys = "PbPb";                                                                  // collision system  (PbPb or pp)
-  int intrate = 50000;                                                                          // interaction rate
+  int intrate = -1;                                                                             // Automatic IR from collision context if -1, else user-defined interaction rate in Hz
   bool flat_gas = true; // if true, the gas density is considered flat in the TPC volume
-  int nFlatGasLoopers = 500;  // number of loopers to be generated per event in case of flat gas
+  unsigned int nFlatGasLoopers = 500;  // number of loopers to be generated per event in case of flat gas
   float fraction_pairs = 0.08; // fraction of loopers
   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

Generators/include/TPCLoopers.h

@@ -91,6 +91,8 @@ class GenTPCLoopers
 
         void SetAdjust(const float &adjust);
 
+        unsigned int getNLoopers() const { return (mNLoopersPairs + mNLoopersCompton); }
+
     private:
         std::unique_ptr<ONNXGenerator> mONNX_pair = nullptr;
         std::unique_ptr<ONNXGenerator> mONNX_compton = nullptr;
@@ -122,6 +124,7 @@ class GenTPCLoopers
         double mTimeLimit = 0.0;                                        // Time limit for the current event
         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
 };
 #endif // GENERATORS_WITH_TPCLOOPERS
 

Generators/src/Generator.cxx

@@ -107,7 +107,7 @@ bool Generator::initLoopersGen()
     }
   }
   const auto& nFlatGasLoopers = loopersParam.nFlatGasLoopers;
-  auto fraction_pairs = loopersParam.fraction_pairs;
+  const auto& fraction_pairs = loopersParam.fraction_pairs;
   std::array<float, 2> multiplier = {loopersParam.multiplier[0], loopersParam.multiplier[1]};
   unsigned int nLoopersPairs = loopersParam.fixedNLoopers[0];
   unsigned int nLoopersCompton = loopersParam.fixedNLoopers[1];
@@ -166,10 +166,6 @@ bool Generator::initLoopersGen()
         LOG(fatal) << "Error: collision system must be either 'PbPb' or 'pp'";
         exit(1);
       } else {
-        if (intrate <= 0) {
-          LOG(fatal) << "Error: interaction rate must be positive!";
-          exit(1);
-        }
         mLoopersGen->SetRate(nclxrate, (colsys == "PbPb") ? true : false, intrate);
         mLoopersGen->SetAdjust(loopersParam.adjust_flatgas);
       }
@@ -217,6 +213,10 @@ Bool_t
       LOG(error) << "Failed to generate loopers event";
       return kFALSE;
     }
+    if (mLoopersGen->getNLoopers() == 0) {
+      LOG(warning) << "No loopers generated for this event";
+      return kTRUE;
+    }
     const auto& looperParticles = mLoopersGen->importParticles();
     if (looperParticles.empty()) {
       LOG(error) << "Failed to import loopers particles";

Generators/src/TPCLoopers.cxx

@@ -437,13 +437,28 @@ void GenTPCLoopers::SetRate(const std::string &rateFile, const bool &isPbPb = tr
     LOG(fatal) << "Error: Could not find fit function '" << fitName << "' in rate file!";
     exit(1);
   }
-  auto ref = static_cast<int>(std::floor(fit->Eval(intRate / 1000.))); // fit expects rate in kHz
+  mInteractionRate = intRate;
+  if (mInteractionRate < 0) {
+    mContextFile = std::filesystem::exists("collisioncontext.root") ? TFile::Open("collisioncontext.root") : nullptr;
+    if (!mContextFile || mContextFile->IsZombie()) {
+      LOG(fatal) << "Error: Interaction rate not provided and collision context file not found!";
+      exit(1);
+    }
+    mCollisionContext = (o2::steer::DigitizationContext*)mContextFile->Get("DigitizationContext");
+    mInteractionRate = std::floor(mCollisionContext->getDigitizerInteractionRate());
+    LOG(info) << "Interaction rate retrieved from collision context: " << mInteractionRate << " Hz";
+    if (mInteractionRate < 0) {
+      LOG(fatal) << "Error: Invalid interaction rate retrieved from collision context!";
+      exit(1);
+    }
+  }
+  auto ref = static_cast<int>(std::floor(fit->Eval(mInteractionRate / 1000.))); // fit expects rate in kHz
   rate_file.Close();
   if (ref <= 0) {
     LOG(fatal) << "Computed flat gas number reference per orbit is <=0";
     exit(1);
   } else {
-    LOG(info) << "Set flat gas number to " << ref << " loopers per orbit using " << fitName << " from " << intRate << " Hz interaction rate.";
+    LOG(info) << "Set flat gas number to " << ref << " loopers per orbit using " << fitName << " from " << mInteractionRate << " Hz interaction rate.";
     auto flat = true;
     setFlatGas(flat, -1, ref);
   }