diff --git a/DataFormats/Parameters/include/DataFormatsParameters/AggregatedRunInfo.h b/DataFormats/Parameters/include/DataFormatsParameters/AggregatedRunInfo.h index bd2cb0c5cbb27..d0347114b5b4c 100644 --- a/DataFormats/Parameters/include/DataFormatsParameters/AggregatedRunInfo.h +++ b/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* ctfFirstRunOrbitVec); + static AggregatedRunInfo buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector* 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); diff --git a/DataFormats/Parameters/src/AggregatedRunInfo.cxx b/DataFormats/Parameters/src/AggregatedRunInfo.cxx index 5495ae73bd6ca..40402a33af68b 100644 --- a/DataFormats/Parameters/src/AggregatedRunInfo.cxx +++ b/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 @@ -42,14 +43,15 @@ o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo_DATA std::map metadata; metadata["runNumber"] = Form("%d", runnumber); auto grpecs = ccdb.getSpecific("GLO/Config/GRPECS", run_mid_timestamp, metadata); + auto grplhcif = ccdb.getSpecific("GLO/Config/GRPLHCIF", run_mid_timestamp); // no run metadata here bool oldFatalState = ccdb.getFatalWhenNull(); ccdb.setFatalWhenNull(false); auto ctp_first_run_orbit = ccdb.getForTimeStamp>("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* ctfFirstRunOrbitVec) +o2::parameters::AggregatedRunInfo AggregatedRunInfo::buildAggregatedRunInfo(int runnumber, long sorMS, long eorMS, long orbitResetMUS, const o2::parameters::GRPECSObject* grpecs, const std::vector* 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 diff --git a/DataFormats/simulation/CMakeLists.txt b/DataFormats/simulation/CMakeLists.txt index fac67cc927562..33c91337c77e9 100644 --- a/DataFormats/simulation/CMakeLists.txt +++ b/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 diff --git a/DataFormats/simulation/include/SimulationDataFormat/InteractionSampler.h b/DataFormats/simulation/include/SimulationDataFormat/InteractionSampler.h index d2ccec147cc4f..47dd4f5e4652d 100644 --- a/DataFormats/simulation/include/SimulationDataFormat/InteractionSampler.h +++ b/DataFormats/simulation/include/SimulationDataFormat/InteractionSampler.h @@ -22,6 +22,7 @@ #include "CommonDataFormat/BunchFilling.h" #include "CommonConstants/LHCConstants.h" #include "MathUtils/RandomRing.h" +#include 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& 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 determineBCIntensityScalesFromHistogram(const TH1F& scales_from_histo); + + const std::vector& getBCIntensityScales() const { return mBCIntensityScales; } + + protected: + int simulateInteractingBC() override; + int getBCJump() const; + + private: + // non-uniformity + std::vector mBCIntensityScales; + ClassDef(NonUniformMuInteractionSampler, 1); +}; + } // namespace steer } // namespace o2 diff --git a/DataFormats/simulation/src/InteractionSampler.cxx b/DataFormats/simulation/src/InteractionSampler.cxx index 5e14e22e5f8db..61b2c4f61bc08 100644 --- a/DataFormats/simulation/src/InteractionSampler.cxx +++ b/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& 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 NonUniformMuInteractionSampler::determineBCIntensityScalesFromHistogram(const TH1F& hist) +{ + std::vector 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; +} \ No newline at end of file diff --git a/DataFormats/simulation/src/SimulationDataLinkDef.h b/DataFormats/simulation/src/SimulationDataLinkDef.h index 15abe9d50390f..8f74bd757e791 100644 --- a/DataFormats/simulation/src/SimulationDataLinkDef.h +++ b/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> + ; diff --git a/DataFormats/simulation/test/testInteractionSampler.cxx b/DataFormats/simulation/test/testInteractionSampler.cxx new file mode 100644 index 0000000000000..b1b3691884ccf --- /dev/null +++ b/DataFormats/simulation/test/testInteractionSampler.cxx @@ -0,0 +1,76 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#define BOOST_TEST_MODULE Test InteractionSampler class +#define BOOST_TEST_MAIN +#define BOOST_TEST_DYN_LINK + +#include +#include "SimulationDataFormat/InteractionSampler.h" +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/AggregatedRunInfo.h" +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "TFile.h" +#include "TGrid.h" +#include + +namespace o2 +{ + +BOOST_AUTO_TEST_CASE(NonUniformSampler) +{ + auto run_number = 559827; + TGrid::Connect("alien"); + if (gGrid) { + auto runInfo = o2::parameters::AggregatedRunInfo::buildAggregatedRunInfo(o2::ccdb::BasicCCDBManager::instance(), run_number); + + o2::steer::NonUniformMuInteractionSampler sampler; + sampler.setBunchFilling(runInfo.grpLHC->getBunchFilling()); + + // the test distribution provided by Igor Altsybeev + auto distr_file = TFile::Open("alien:///alice/cern.ch/user/s/swenzel/AliceO2_TestData/NBcVTX_559827/hBcTVX_data_PbPb_24ar_559827.root"); + + // + if (distr_file && !distr_file->IsZombie()) { + auto hist = distr_file->Get("hBcTVX"); + if (hist) { + sampler.init(); + sampler.setBCIntensityScales(*hist); + + // sample into a vector of a certain size + std::vector samples; + + int N = 100000; + samples.resize(N); + + sampler.generateCollisionTimes(samples); + + // fill an output histogram + auto output_hist = (TH1F*)hist->Clone("h2"); // make a full copy + output_hist->Reset(); + + for (const auto& sample : samples) { + output_hist->Fill(sample.bc); + } + + // Write out + auto fout = TFile::Open("NBCVTX_out.root", "RECREATE"); + fout->WriteObject(output_hist, "NBcVTX"); + fout->Close(); + + // compare mean values of original and newly sampled hist + BOOST_CHECK_CLOSE(hist->GetMean(), output_hist->GetMean(), 0.5); + } + } + } +} + +} // namespace o2