Merge 6e4b06d into sapling-pr-archive-ktf

Author: ktf

DataFormats/Detectors/TRD/include/DataFormatsTRD/Digit.h

@@ -59,7 +59,7 @@ class Digit
   Digit(int det, int row, int pad, ArrayADC adc, int phase = 0);
   Digit(int det, int row, int pad); // add adc data and pretrigger phase in a separate step
   Digit(int det, int rob, int mcm, int channel, ArrayADC adc, int phase = 0);
-  Digit(int det, int rob, int mcm, int channel); // add adc data in a seperate step
+  Digit(int det, int rob, int mcm, int channel, int phase = 0); // add adc data
 
   // Copy
   Digit(const Digit&) = default;
@@ -74,9 +74,11 @@ class Digit
   void setDetector(int det) { mDetector = ((mDetector & 0xf000) | (det & 0xfff)); }
   void setADC(ArrayADC const& adc) { mADC = adc; }
   void setADC(const gsl::span<ADC_t>& adc) { std::copy(adc.begin(), adc.end(), mADC.begin()); }
-  void setPreTrigPhase(int phase) { mDetector = (((phase & 0xf) << 12) | (mDetector & 0xfff)); }
+  // set the trigger phase make sure it is mapped to 2 bits as it can only have 4 valid numbers shifted 0,3,6,9 or 1,4,7,10 etc.
+  void setPreTrigPhase(int phase);
   // Get methods
   int getDetector() const { return mDetector & 0xfff; }
+  int getDetectorInFull() const { return mDetector; } // return the entire mDetector 16 bits, so far only for CTF encoding.
   int getHCId() const { return (mDetector & 0xfff) * 2 + (mROB % 2); }
   int getPadRow() const { return HelperMethods::getPadRowFromMCM(mROB, mMCM); }
   int getPadCol() const { return HelperMethods::getPadColFromADC(mROB, mMCM, mChannel); }

DataFormats/Detectors/TRD/include/DataFormatsTRD/PHData.h

@@ -61,6 +61,66 @@ class PHData
 
   ClassDefNV(PHData, 1);
 };
+
+/*
+  This data type is used to send around the information required to fill PH plots per chamber
+
+  |19|18|17|16|15|14|13|12|11|10|09|08|07|06|05|04|03|02|01|00|
+  -------------------------------------------------------------
+  |type |nNeighb |   time bin   |        detector number      |
+  -------------------------------------------------------------
+*/
+/*
+  This data type is used to send around the information required to fill PH plots per chamber
+
+  |15|14|13|12|11|10|09|08|07|06|05|04|03|02|01|00|
+  ------------------------------------------------
+  |                 ADC sum for all neigbours     |
+  ------------------------------------------------
+*/
+
+class PHDataHD
+{
+ public:
+  enum Origin : uint8_t {
+    ITSTPCTRD,
+    TPCTRD,
+    TRACKLET,
+    OTHER
+  };
+
+  PHDataHD() = default;
+  PHDataHD(int adc, int det, int tb, int nb, int type) { set(adc, det, tb, nb, type); }
+
+  void set(int adc, int det, int tb, int nb, int type)
+  {
+    mDetector = det;
+    mTimeBin = tb;
+    mType = type;
+    mNNeighbours = nb;
+    mADC = adc;
+  }
+
+  // the ADC sum for given time bin for up to three neighbours
+  int getADC() const { return mADC; }
+  // the TRD detector number
+  int getDetector() const { return mDetector; }
+  // the given time bin
+  int getTimebin() const { return mTimeBin; }
+  // number of neighbouring digits for which the ADC is accumulated
+  int getNNeighbours() const { return mNNeighbours; }
+  // the origin of this point: digit on ITS-TPC-TRD track, ... (see enum Origin above)
+  int getType() const { return mType; }
+
+ private:
+  uint16_t mDetector{0};
+  uint8_t mTimeBin{0};
+  uint8_t mType{0};
+  uint8_t mNNeighbours{0};
+  uint16_t mADC{0};
+
+  ClassDefNV(PHDataHD, 1);
+};
 } // namespace o2::trd
 
 #endif // ALICEO2_TRD_PHDATA_H_

DataFormats/Detectors/TRD/src/DataFormatsTRDLinkDef.h

@@ -43,6 +43,7 @@
 #pragma link C++ class o2::trd::ChannelInfo + ;
 #pragma link C++ class o2::trd::ChannelInfoContainer + ;
 #pragma link C++ struct o2::trd::PHData + ;
+#pragma link C++ struct o2::trd::PHDataHD + ;
 #pragma link C++ class o2::trd::TRDDataCountersPerTimeFrame + ;
 #pragma link C++ class o2::trd::DataCountersPerTrigger + ;
 #pragma link C++ class std::vector < o2::trd::Tracklet64> + ;
@@ -56,6 +57,7 @@
 #pragma link C++ class std::vector < o2::trd::GainCalibHistos> + ;
 #pragma link C++ class std::vector < o2::trd::T0FitHistos> + ;
 #pragma link C++ class std::vector < o2::trd::PHData> + ;
+#pragma link C++ class std::vector < o2::trd::PHDataHD> + ;
 #pragma link C++ class std::vector < o2::trd::KrCluster> + ;
 #pragma link C++ class std::vector < o2::trd::KrClusterTriggerRecord> + ;
 #pragma link C++ class std::vector < o2::trd::DataCountersPerTrigger> + ;

DataFormats/Detectors/TRD/src/Digit.cxx

@@ -12,6 +12,7 @@
 #include "DataFormatsTRD/Digit.h"
 #include <iostream>
 #include <algorithm>
+#include "fairlogger/Logger.h"
 
 namespace o2::trd
 {
@@ -46,12 +47,18 @@ Digit::Digit(int det, int rob, int mcm, int channel, ArrayADC adc, int pretrigph
   setPreTrigPhase(pretrigphase);
 }
 
-Digit::Digit(int det, int rob, int mcm, int channel) // add adc data in a seperate step
+Digit::Digit(int det, int rob, int mcm, int channel, int pretrigphase) // add adc data in a seperate step
 {
   setDetector(det);
   setROB(rob);
   setMCM(mcm);
   setChannel(channel);
+  setPreTrigPhase(pretrigphase);
+}
+
+void Digit::setPreTrigPhase(int phase)
+{
+  mDetector = ((((phase) & 0x3) << 12) | (mDetector & 0xfff));
 }
 
 bool Digit::isSharedDigit() const

DataFormats/Parameters/include/DataFormatsParameters/AggregatedRunInfo.h

@@ -23,6 +23,7 @@ namespace o2::parameters
 {
 
 class GRPECSObject;
+class GRPLHCIFData;
 
 /// Composite struct where one may collect important global properties of data "runs"
 /// aggregated from various sources (GRPECS, RunInformation CCDB entries, etc.).
@@ -39,8 +40,9 @@ struct AggregatedRunInfo {
 
   // we may have pointers to actual data source objects GRPECS, ...
   const o2::parameters::GRPECSObject* grpECS = nullptr; // pointer to GRPECSobject (fetched during struct building)
+  const o2::parameters::GRPLHCIFData* grpLHC = nullptr;
 
-  static AggregatedRunInfo buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector<Long64_t>* ctfFirstRunOrbitVec);
+  static AggregatedRunInfo buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector<Long64_t>* ctfFirstRunOrbitVec, const o2::parameters::GRPLHCIFData* grplhcif = nullptr);
 
   // fills and returns AggregatedRunInfo for a given data run number.
   static AggregatedRunInfo buildAggregatedRunInfo_DATA(o2::ccdb::CCDBManagerInstance& ccdb, int runnumber);

DataFormats/Parameters/src/AggregatedRunInfo.cxx

@@ -15,6 +15,7 @@
 #include "DataFormatsParameters/AggregatedRunInfo.h"
 #include "CCDB/BasicCCDBManager.h"
 #include "DataFormatsParameters/GRPECSObject.h"
+#include "DataFormatsParameters/GRPLHCIFData.h"
 #include "CommonConstants/LHCConstants.h"
 #include "Framework/Logger.h"
 #include <map>
@@ -42,14 +43,15 @@ o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo_DATA
   std::map<std::string, std::string> metadata;
   metadata["runNumber"] = Form("%d", runnumber);
   auto grpecs = ccdb.getSpecific<o2::parameters::GRPECSObject>("GLO/Config/GRPECS", run_mid_timestamp, metadata);
+  auto grplhcif = ccdb.getSpecific<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", run_mid_timestamp); // no run metadata here
   bool oldFatalState = ccdb.getFatalWhenNull();
   ccdb.setFatalWhenNull(false);
   auto ctp_first_run_orbit = ccdb.getForTimeStamp<std::vector<Long64_t>>("CTP/Calib/FirstRunOrbit", run_mid_timestamp);
   ccdb.setFatalWhenNull(oldFatalState);
-  return buildAggregatedRunInfo(runnumber, sor, eor, tsOrbitReset, grpecs, ctp_first_run_orbit);
+  return buildAggregatedRunInfo(runnumber, sor, eor, tsOrbitReset, grpecs, ctp_first_run_orbit, grplhcif);
 }
 
-o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector<Long64_t>* ctfFirstRunOrbitVec)
+o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector<Long64_t>* ctfFirstRunOrbitVec, const o2::parameters::GRPLHCIFData* grplhcif)
 {
   auto nOrbitsPerTF = grpecs->getNHBFPerTF();
   // calculate SOR/EOR orbits
@@ -81,7 +83,7 @@ o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo(int
       orbitSOR = (orbitSOR / nOrbitsPerTF + 1) * nOrbitsPerTF;
     }
   }
-  return AggregatedRunInfo{runnumber, sorMS, eorMS, nOrbitsPerTF, orbitResetMUS, orbitSOR, orbitEOR, grpecs};
+  return AggregatedRunInfo{runnumber, sorMS, eorMS, nOrbitsPerTF, orbitResetMUS, orbitSOR, orbitEOR, grpecs, grplhcif};
 }
 
 namespace

DataFormats/simulation/CMakeLists.txt

@@ -55,6 +55,11 @@ o2_target_root_dictionary(
 # * src/SimulationDataLinkDef.h
 # * and not src/SimulationDataFormatLinkDef.h
 
+o2_add_test(InteractionSampler
+            SOURCES test/testInteractionSampler.cxx
+            COMPONENT_NAME SimulationDataFormat
+            PUBLIC_LINK_LIBRARIES O2::SimulationDataFormat)
+
 o2_add_test(BasicHits
             SOURCES test/testBasicHits.cxx
             COMPONENT_NAME SimulationDataFormat

DataFormats/simulation/include/SimulationDataFormat/InteractionSampler.h

@@ -22,6 +22,7 @@
 #include "CommonDataFormat/BunchFilling.h"
 #include "CommonConstants/LHCConstants.h"
 #include "MathUtils/RandomRing.h"
+#include <TH1F.h>
 
 namespace o2
 {
@@ -130,6 +131,30 @@ class FixedSkipBC_InteractionSampler : public InteractionSampler
   ClassDef(FixedSkipBC_InteractionSampler, 1);
 };
 
+// A version of the interaction sampler which can sample according to non-uniform mu(bc) as
+// observed during data taking.
+class NonUniformMuInteractionSampler : public InteractionSampler
+{
+ public:
+  NonUniformMuInteractionSampler() : InteractionSampler() { mBCIntensityScales.resize(o2::constants::lhc::LHCMaxBunches, 1); }
+  bool setBCIntensityScales(const std::vector<float>& scales_from_vector);
+  bool setBCIntensityScales(const TH1F& scales_from_histo); // initialize scales
+
+  // helper function to determine the scales from a histogram (count from event selection analysis)
+  std::vector<float> determineBCIntensityScalesFromHistogram(const TH1F& scales_from_histo);
+
+  const std::vector<float>& getBCIntensityScales() const { return mBCIntensityScales; }
+
+ protected:
+  int simulateInteractingBC() override;
+  int getBCJump() const;
+
+ private:
+  // non-uniformity
+  std::vector<float> mBCIntensityScales;
+  ClassDef(NonUniformMuInteractionSampler, 1);
+};
+
 } // namespace steer
 } // namespace o2
 

DataFormats/simulation/src/InteractionSampler.cxx

@@ -115,8 +115,8 @@ const o2::InteractionTimeRecord& InteractionSampler::generateCollisionTime()
 int InteractionSampler::simulateInteractingBC()
 {
   // Returns number of collisions assigned to selected BC
-
   nextCollidingBC(mBCJumpGenerator.getNextValue());
+
   // once BC is decided, enforce at least one interaction
   int ncoll = mNCollBCGenerator.getNextValue();
 
@@ -162,3 +162,98 @@ void InteractionSampler::setBunchFilling(const std::string& bcFillingFile)
   mBCFilling = *bc;
   delete bc;
 }
+
+// ________________________________________________
+bool NonUniformMuInteractionSampler::setBCIntensityScales(const std::vector<float>& scales_from_vector)
+{
+  // Sets the intensity scales per bunch crossing index
+  // The length of this vector needs to be compatible with the bunch filling chosen
+  mBCIntensityScales = scales_from_vector;
+
+  if (scales_from_vector.size() != mInteractingBCs.size()) {
+    LOG(error) << "Scaling factors and bunch filling scheme are not compatible. Not doing anything";
+    return false;
+  }
+
+  float sum = 0.;
+  for (auto v : mBCIntensityScales) {
+    sum += std::abs(v);
+  }
+  if (sum == 0) {
+    LOGP(warn, "total intensity is 0, assuming uniform");
+    for (auto& v : mBCIntensityScales) {
+      v = 1.f;
+    }
+  } else { // normalize
+    float norm = mBCIntensityScales.size() / sum;
+    for (auto& v : mBCIntensityScales) {
+      v = std::abs(v) * norm;
+    }
+  }
+  return false;
+}
+
+// ________________________________________________
+
+bool NonUniformMuInteractionSampler::setBCIntensityScales(const TH1F& hist)
+{
+  return setBCIntensityScales(determineBCIntensityScalesFromHistogram(hist));
+}
+
+std::vector<float> NonUniformMuInteractionSampler::determineBCIntensityScalesFromHistogram(const TH1F& hist)
+{
+  std::vector<float> scales;
+  // we go through the BCs and query the count from histogram
+  for (auto bc : mInteractingBCs) {
+    scales.push_back(hist.GetBinContent(bc + 1));
+  }
+  return scales;
+}
+
+int NonUniformMuInteractionSampler::getBCJump() const
+{
+  auto muFunc = [this](int bc_position) {
+    return mBCIntensityScales[bc_position % mInteractingBCs.size()] * mMuBC;
+  };
+
+  double U = gRandom->Rndm();    // uniform (0,1)
+  double T = -std::log(1.0 - U); // threshold
+  double sumMu = 0.0;
+  int offset = 0;
+  auto bcStart = mCurrBCIdx; // the current bc
+
+  while (sumMu < T) {
+    auto mu_here = muFunc(bcStart + offset); // mu at next BC
+    sumMu += mu_here;
+    if (sumMu >= T) {
+      break; // found BC with at least one collision
+    }
+    ++offset;
+  }
+  return offset;
+}
+
+int NonUniformMuInteractionSampler::simulateInteractingBC()
+{
+  nextCollidingBC(getBCJump());
+
+  auto muFunc = [this](int bc_position) {
+    return mBCIntensityScales[bc_position % mInteractingBCs.size()] * mMuBC;
+  };
+
+  // now sample number of collisions in chosenBC, conditioned >=1:
+  double mu_chosen = muFunc(mCurrBCIdx); // or does it need to be mCurrBCIdx
+  int ncoll = 0;
+  do {
+    ncoll = gRandom->Poisson(mu_chosen);
+  } while (ncoll == 0);
+
+  // assign random time withing a bunch
+  for (int i = ncoll; i--;) {
+    mTimeInBC.push_back(mCollTimeGenerator.getNextValue());
+  }
+  if (ncoll > 1) { // sort in DECREASING time order (we are reading vector from the end)
+    std::sort(mTimeInBC.begin(), mTimeInBC.end(), [](const float a, const float b) { return a > b; });
+  }
+  return ncoll;
+}

DataFormats/simulation/src/SimulationDataLinkDef.h

@@ -25,6 +25,7 @@
 
 #pragma link C++ class o2::steer::InteractionSampler + ;
 #pragma link C++ class o2::steer::FixedSkipBC_InteractionSampler + ;
+#pragma link C++ class o2::steer::NonUniformMuInteractionSampler + ;
 #pragma link C++ class o2::sim::StackParam + ;
 #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::sim::StackParam> + ;
 #pragma link C++ class o2::MCTrackT < double> + ;