From 72ea0a9833d12bc07e14c2a1382213f9fff7921b Mon Sep 17 00:00:00 2001
From: Volker Friese <v.friese@gsi.de>
Date: Mon, 26 Jun 2023 12:42:38 +0000
Subject: [PATCH] Introduce possibility to simulate time-based with constant
time between events. Refs #1792.
---
macro/mcbm/mcbm_digi.C | 4 +-
macro/mvd/qa/mvd_qa3_digitize.C | 2 +-
macro/run/run_digi.C | 8 ++--
sim/response/CMakeLists.txt | 3 ++
sim/response/base/CbmDigitization.cxx | 38 ++++++------------
sim/response/base/CbmDigitization.h | 44 +++++++++++----------
sim/response/base/CbmDigitizationConfig.cxx | 4 +-
sim/response/base/CbmDigitizationSource.cxx | 15 ++++---
sim/response/base/CbmDigitizationSource.h | 7 ++--
sim/response/base/CbmMCInputSet.cxx | 30 ++++++++++----
sim/response/base/CbmMCInputSet.h | 17 ++++----
sim/response/base/Defs.h | 37 +++++++++++++++++
12 files changed, 131 insertions(+), 78 deletions(-)
create mode 100644 sim/response/base/Defs.h
diff --git a/macro/mcbm/mcbm_digi.C b/macro/mcbm/mcbm_digi.C
index 3cc93000e9..1ac53fac74 100644
--- a/macro/mcbm/mcbm_digi.C
+++ b/macro/mcbm/mcbm_digi.C
@@ -68,11 +68,11 @@ void mcbm_digi(Int_t nEvents = 10, // Number of events to p
// ----- Digitization run ---------------------------------------------
CbmDigitization run;
- run.AddInput(inFile, eventRate);
+ run.AddInput(inFile, cbm::sim::TimeDist::Poisson, eventRate);
run.SetOutputFile(outFile, overwrite);
run.SetParameterRootFile(parFile);
run.SetTimeSliceLength(timeSliceLength);
- run.SetEventMode(eventMode);
+ run.SetMode(eventMode ? cbm::sim::Mode::EventByEvent : cbm::sim::Mode::Timebased);
run.Run(nEvents);
// ------------------------------------------------------------------------
diff --git a/macro/mvd/qa/mvd_qa3_digitize.C b/macro/mvd/qa/mvd_qa3_digitize.C
index 6cb7ef0349..d2d8816c9c 100644
--- a/macro/mvd/qa/mvd_qa3_digitize.C
+++ b/macro/mvd/qa/mvd_qa3_digitize.C
@@ -99,7 +99,7 @@ void mvd_qa3_digitize(const char* setup = "sis100_electron")
run.AddInput(inFile);
run.SetOutputFile(outFile, kTRUE);
run.SetParameterRootFile(parFile);
- run.SetEventMode();
+ run.SetMode(cbm::sim::Mode::EventByEvent);
run.SetDigitizer(ECbmModuleId::kMvd, digi);
diff --git a/macro/run/run_digi.C b/macro/run/run_digi.C
index bede1dc98d..880d1b0f46 100644
--- a/macro/run/run_digi.C
+++ b/macro/run/run_digi.C
@@ -129,16 +129,18 @@ void run_digi(TString inputEvents = "", Int_t nEvents = -1, TString output = "",
// ----- Digitization run ---------------------------------------------
CbmDigitization run;
- run.AddInput(0, evntFile, eventRate);
+ cbm::sim::Mode mode = (eventRate < 0. ? cbm::sim::Mode::EventByEvent : cbm::sim::Mode::Timebased);
+ cbm::sim::TimeDist timeDist = cbm::sim::TimeDist::Uniform;
+ run.AddInput(0, evntFile, timeDist, eventRate);
if (!eventMode) {
if (!inputSignal.IsNull()) run.EmbedInput(1, signFile, 0);
- if (!inputBeam.IsNull()) run.AddInput(2, beamFile, beamRate);
+ if (!inputBeam.IsNull()) run.AddInput(2, beamFile, timeDist, beamRate);
}
run.SetOutputFile(outFile, overwrite);
run.SetMonitorFile(moniFile);
run.SetParameterRootFile(parFile);
run.SetTimeSliceLength(tsLength);
- run.SetEventMode(eventMode);
+ run.SetMode(mode);
run.SetProduceNoise(kFALSE);
if (nEvents < 0) run.Run();
else
diff --git a/sim/response/CMakeLists.txt b/sim/response/CMakeLists.txt
index dde4f0799b..88992d439d 100644
--- a/sim/response/CMakeLists.txt
+++ b/sim/response/CMakeLists.txt
@@ -40,3 +40,6 @@ set(INTERFACE_DEPENDENCIES
)
generate_cbm_library()
+
+install(FILES base/Defs.h DESTINATION include/)
+
diff --git a/sim/response/base/CbmDigitization.cxx b/sim/response/base/CbmDigitization.cxx
index 3acf0c1c6e..9080ba6870 100644
--- a/sim/response/base/CbmDigitization.cxx
+++ b/sim/response/base/CbmDigitization.cxx
@@ -37,28 +37,12 @@
#include <cassert>
+using cbm::sim::Mode;
+using cbm::sim::TimeDist;
+
// ----- Constructor ----------------------------------------------------
-CbmDigitization::CbmDigitization()
- : TNamed("CbmDigitization", "Digitisation Run")
- , fIsInit(kFALSE)
- , fEventMode(kFALSE)
- , fTimeSliceLength(-1.)
- , fProduceNoise(kTRUE)
- , fCreateMatches(kTRUE)
- , fDigitizers()
- , fDaq(new CbmDaq())
- , fSource(new CbmDigitizationSource())
- , fOutFile()
- , fParRootFile()
- , fParAsciiFiles()
- , fMoniFile()
- , fOverwriteOutput(kFALSE)
- , fGenerateRunInfo(kTRUE)
- , fRun(0)
-{
- SetDefaultBranches();
-}
+CbmDigitization::CbmDigitization() : TNamed("CbmDigitization", "Digitisation Run") { SetDefaultBranches(); }
// --------------------------------------------------------------------------
@@ -75,14 +59,14 @@ CbmDigitization::~CbmDigitization()
// ----- Add an input file ----------------------------------------------
-void CbmDigitization::AddInput(UInt_t inputId, TString fileName, Double_t eventRate, ECbmTreeAccess mode)
+void CbmDigitization::AddInput(UInt_t inputId, TString fileName, TimeDist dist, Double_t eventRate, ECbmTreeAccess mode)
{
if (gSystem->AccessPathName(fileName)) LOG(fatal) << fName << ": input file " << fileName << " does not exist!";
if (mode != ECbmTreeAccess::kRegular)
LOG(fatal) << fName << ": access modes other than kRegular are not yet supported!";
TChain* chain = new TChain("cbmsim");
chain->Add(fileName.Data());
- fSource->AddInput(inputId, chain, eventRate, mode);
+ fSource->AddInput(inputId, chain, dist, eventRate, mode);
}
// --------------------------------------------------------------------------
@@ -160,7 +144,7 @@ Int_t CbmDigitization::CreateDefaultDigitizers()
if (it->second->GetDigitizer() != nullptr) continue;
// --- Skip MVD for time-based mode
- if (it->first == ECbmModuleId::kMvd && (!fEventMode)) {
+ if (it->first == ECbmModuleId::kMvd && fMode == Mode::Timebased) {
LOG(info) << "MVD digitizer is not available "
<< "in time-based mode. ";
continue;
@@ -362,13 +346,13 @@ void CbmDigitization::Run(Int_t event1, Int_t event2)
if (fGenerateRunInfo) LOG(info) << fName << ": Run info will be generated.";
// --- Create DAQ
- if (fEventMode) fDaq = new CbmDaq(kTRUE);
+ if (fMode == Mode::EventByEvent) fDaq = new CbmDaq(kTRUE);
else
fDaq = new CbmDaq(fTimeSliceLength);
// --- Register source
- if (fEventMode) fSource->SetEventMode();
+ fSource->SetMode(fMode);
run->SetSource(fSource);
@@ -391,7 +375,7 @@ void CbmDigitization::Run(Int_t event1, Int_t event2)
CbmDigitizeBase* digitizer = it->second->GetDigitizer();
if (it->second->IsActive() && digitizer != nullptr) {
fDaq->SetDigitizer(it->first, digitizer);
- if (fEventMode) digitizer->SetEventMode();
+ if (fMode == Mode::EventByEvent) digitizer->SetEventMode();
digitizer->SetProduceNoise(fProduceNoise);
digitizer->SetCreateMatches(fCreateMatches);
run->AddTask(digitizer);
@@ -401,7 +385,7 @@ void CbmDigitization::Run(Int_t event1, Int_t event2)
// --- In event-by-event mode: also empty events (time-slices) are stored
- if (fEventMode) StoreAllTimeSlices();
+ if (fMode == Mode::EventByEvent) StoreAllTimeSlices();
// --- Register DAQ
diff --git a/sim/response/base/CbmDigitization.h b/sim/response/base/CbmDigitization.h
index bea69df204..e2c3cfcd0d 100644
--- a/sim/response/base/CbmDigitization.h
+++ b/sim/response/base/CbmDigitization.h
@@ -25,6 +25,8 @@
#include <map>
#include <vector>
+#include "Defs.h"
+
class TGeoManager;
class CbmDigitizeBase;
@@ -46,8 +48,8 @@ public:
** @param eventRate Rate for events from input file [1/s]
** @param mode Tree access mode (kRegular / kRepeat / kRandom)
**/
- void AddInput(UInt_t inputId, TString fileName, Double_t eventRate = -1.,
- ECbmTreeAccess mode = ECbmTreeAccess::kRegular);
+ void AddInput(UInt_t inputId, TString fileName, cbm::sim::TimeDist dist = cbm::sim::TimeDist::Poisson,
+ Double_t eventRate = 0., ECbmTreeAccess mode = ECbmTreeAccess::kRegular);
/** @brief Add an input file
@@ -58,9 +60,10 @@ public:
** Shortcut for legacy reasons, when only one input file is used.
** This will set the inputId to zero. Repeated use will lead to abort.
**/
- void AddInput(TString fileName, Double_t eventRate = -1., ECbmTreeAccess mode = ECbmTreeAccess::kRegular)
+ void AddInput(TString fileName, cbm::sim::TimeDist dist = cbm::sim::TimeDist::Poisson, Double_t eventRate = 0.,
+ ECbmTreeAccess mode = ECbmTreeAccess::kRegular)
{
- AddInput(0, fileName, eventRate, mode);
+ AddInput(0, fileName, dist, eventRate, mode);
}
/** @brief Add an ASCII parameter file
@@ -156,7 +159,7 @@ public:
** In the event-by-event mode, one time slice will be created for
** each input event. There will be no interference between events.
**/
- void SetEventMode(Bool_t choice = kTRUE) { fEventMode = choice; }
+ void SetMode(cbm::sim::Mode mode) { fMode = mode; }
/** @brief Set the monitor file name
@@ -230,21 +233,22 @@ public:
private:
- Bool_t fIsInit;
- Bool_t fEventMode;
- Double_t fTimeSliceLength;
- Bool_t fProduceNoise;
- Bool_t fCreateMatches;
- std::map<ECbmModuleId, CbmDigitizeInfo*> fDigitizers; //!
- CbmDaq* fDaq;
- CbmDigitizationSource* fSource; ///< Input source
- TString fOutFile; ///< Output data (digis)
- TString fParRootFile; ///< ROOT parameter file
- TList fParAsciiFiles; ///< ASCII parameter files
- TString fMoniFile; ///< Resource monitoring information
- Bool_t fOverwriteOutput;
- Bool_t fGenerateRunInfo;
- Int_t fRun;
+ Bool_t fIsInit = kFALSE;
+ cbm::sim::Mode fMode = cbm::sim::Mode::Timebased;
+ Double_t fTimeSliceLength = -1.;
+ Bool_t fProduceNoise = kTRUE;
+ ;
+ Bool_t fCreateMatches = kTRUE;
+ std::map<ECbmModuleId, CbmDigitizeInfo*> fDigitizers = {}; //!
+ CbmDaq* fDaq = new CbmDaq();
+ CbmDigitizationSource* fSource = new CbmDigitizationSource(); ///< Input source
+ TString fOutFile = {}; ///< Output data (digis)
+ TString fParRootFile = {}; ///< ROOT parameter file
+ TList fParAsciiFiles = {}; ///< ASCII parameter files
+ TString fMoniFile = {}; ///< Resource monitoring information
+ Bool_t fOverwriteOutput = kFALSE;
+ Bool_t fGenerateRunInfo = kTRUE;
+ Int_t fRun = 0;
/** @brief Copy constructor forbidden **/
diff --git a/sim/response/base/CbmDigitizationConfig.cxx b/sim/response/base/CbmDigitizationConfig.cxx
index 980771d232..dd4422c4f2 100644
--- a/sim/response/base/CbmDigitizationConfig.cxx
+++ b/sim/response/base/CbmDigitizationConfig.cxx
@@ -90,7 +90,7 @@ bool CbmDigitizationConfig::SetIO(CbmDigitization& obj, const pt::ptree& moduleT
}
float rate = pt_input.get<float>("rate", -1.);
LOG(info) << "CbmDigitizationConfig: Adding input: " << traFileName;
- obj.AddInput(id, path + ".tra.root", rate, treeAccessMode);
+ obj.AddInput(id, path + ".tra.root", cbm::sim::TimeDist::Poisson, rate, treeAccessMode);
}
if (parameterSource) {
if (!parametersSet) {
@@ -134,7 +134,7 @@ bool CbmDigitizationConfig::SetDigitizationParameters(CbmDigitization& obj, cons
return false;
}
}
- obj.SetEventMode(eventMode);
+ obj.SetMode(eventMode ? cbm::sim::Mode::EventByEvent : cbm::sim::Mode::Timebased);
if (timeSliceLength) obj.SetTimeSliceLength(timeSliceLength.get());
if (startTime) obj.SetStartTime(startTime.get());
diff --git a/sim/response/base/CbmDigitizationSource.cxx b/sim/response/base/CbmDigitizationSource.cxx
index 362a59e434..ca060b36b0 100644
--- a/sim/response/base/CbmDigitizationSource.cxx
+++ b/sim/response/base/CbmDigitizationSource.cxx
@@ -22,6 +22,9 @@
#include <iomanip>
#include <utility>
+using cbm::sim::Mode;
+using cbm::sim::TimeDist;
+
// ----- Constructor -----------------------------------------------------
CbmDigitizationSource::CbmDigitizationSource()
@@ -38,7 +41,7 @@ CbmDigitizationSource::CbmDigitizationSource()
, fCurrentInputId(0)
, fCurrentRunId(0)
, fFirstCall(kTRUE)
- , fEventMode(kFALSE)
+ , fMode(cbm::sim::Mode::Timebased)
, fCurrentInputSet(nullptr)
, fSwitchInputSet(kFALSE)
{
@@ -70,7 +73,7 @@ Bool_t CbmDigitizationSource::ActivateObject(TObject** object, const char* branc
// ----- Add a transport input -------------------------------------------
-void CbmDigitizationSource::AddInput(UInt_t inputId, TChain* chain, Double_t rate, ECbmTreeAccess mode)
+void CbmDigitizationSource::AddInput(UInt_t inputId, TChain* chain, TimeDist dist, Double_t rate, ECbmTreeAccess mode)
{
// Catch input ID already being used
@@ -80,7 +83,7 @@ void CbmDigitizationSource::AddInput(UInt_t inputId, TChain* chain, Double_t rat
} //? Input ID already in use
// Create a new inputSet and add the input to it
- CbmMCInputSet* inputSet = new CbmMCInputSet(rate);
+ CbmMCInputSet* inputSet = new CbmMCInputSet(dist, rate);
inputSet->AddInput(inputId, chain, mode);
// If it is the first input set, it defines the reference branch list.
@@ -217,7 +220,7 @@ Bool_t CbmDigitizationSource::Init()
LOG(fatal) << "DigitizationSource: No input sets defined!";
return kFALSE;
}
- if (fEventMode && fInputMap.size() != 1) {
+ if (fMode == Mode::EventByEvent && fInputMap.size() != 1) {
LOG(fatal) << "DigitizationSource: More than one input defined "
<< "in event-by-event mode!";
return kFALSE;
@@ -246,7 +249,7 @@ Bool_t CbmDigitizationSource::Init()
ActivateObject(object, "MCEventHeader.");
// Set the time of the first event for each input set
- if (!fEventMode) {
+ if (fMode == Mode::Timebased) {
for (size_t iSet = 0; iSet < fInputSets.size(); iSet++) {
CbmMCInputSet* inputSet = fInputSets.at(iSet);
Double_t time = fTimeStart + inputSet->GetDeltaT();
@@ -282,7 +285,7 @@ Int_t CbmDigitizationSource::ReadEvent(UInt_t event)
// In the event-by-event mode, get the respective event from the first
// input; the event time is fStartTime.
- if (fEventMode) return ReadEventByEvent(event);
+ if (fMode == Mode::EventByEvent) return ReadEventByEvent(event);
// If the last used input set was exhausted, switch to a the next one
if (fSwitchInputSet) {
diff --git a/sim/response/base/CbmDigitizationSource.h b/sim/response/base/CbmDigitizationSource.h
index 2a4d3f20e6..17517920ce 100644
--- a/sim/response/base/CbmDigitizationSource.h
+++ b/sim/response/base/CbmDigitizationSource.h
@@ -75,7 +75,8 @@ public:
** @param chain Pointer to input chain
** @param mode Tree access mode (kRegular / kRepeat / kRandom)
**/
- void AddInput(UInt_t inputId, TChain* chain, Double_t rate, ECbmTreeAccess mode = ECbmTreeAccess::kRegular);
+ void AddInput(UInt_t inputId, TChain* chain, cbm::sim::TimeDist dist, Double_t rate,
+ ECbmTreeAccess mode = ECbmTreeAccess::kRegular);
/** @brief Maximal entry number the source can run to
@@ -190,7 +191,7 @@ public:
** In the event-by-event mode, only the first input is processed.
** No event start time is generated; the event time is always zero.
**/
- void SetEventMode(Bool_t choice = kTRUE) { fEventMode = choice; }
+ void SetMode(cbm::sim::Mode mode) { fMode = mode; }
/** @brief Abstract in base class. No implementation here.
@@ -226,7 +227,7 @@ private:
Int_t fCurrentInputId;
Int_t fCurrentRunId;
Bool_t fFirstCall;
- Bool_t fEventMode;
+ cbm::sim::Mode fMode;
CbmMCInputSet* fCurrentInputSet;
Bool_t fSwitchInputSet; // Flag to switch the input set at next ReadEvent
diff --git a/sim/response/base/CbmMCInputSet.cxx b/sim/response/base/CbmMCInputSet.cxx
index 3cc0d400ff..74a4be2b68 100644
--- a/sim/response/base/CbmMCInputSet.cxx
+++ b/sim/response/base/CbmMCInputSet.cxx
@@ -14,24 +14,27 @@
#include <cassert>
+using cbm::sim::TimeDist;
+
// ----- Default constructor ---------------------------------------------
-CbmMCInputSet::CbmMCInputSet() : CbmMCInputSet(-1.) {}
+CbmMCInputSet::CbmMCInputSet() : CbmMCInputSet(TimeDist::Poisson, -1.) {}
// ---------------------------------------------------------------------------
// ----- Constructor -----------------------------------------------------
-CbmMCInputSet::CbmMCInputSet(Double_t rate)
+CbmMCInputSet::CbmMCInputSet(TimeDist dist, Double_t eventRate)
: TObject()
- , fRate(rate)
+ , fTimeDist(dist)
+ , fEventRate(eventRate)
, fInputs()
, fInputHandle()
, fBranches()
, fDeltaDist(nullptr)
{
- if (rate > 0.) {
- Double_t mean = 1.e9 / rate; // mean time between events
+ if (fTimeDist == TimeDist::Poisson && eventRate > 0.) {
+ Double_t mean = 1.e9 / eventRate; // mean time between events in ns
fDeltaDist = new TF1("DeltaDist", "exp(-x/[0])/[0]", 0., 10. * mean, "NL");
fDeltaDist->SetParameter(0, mean);
}
@@ -140,8 +143,21 @@ Bool_t CbmMCInputSet::CheckBranchList(CbmMCInput* input)
// ----- Time difference to next event -----------------------------------
Double_t CbmMCInputSet::GetDeltaT()
{
- if (!fDeltaDist) return 0.;
- return fDeltaDist->GetRandom();
+
+ // --- Poisson distribution: delta_t sampled from exponential distribution
+ if (fTimeDist == TimeDist::Poisson) {
+ assert(fDeltaDist);
+ return fDeltaDist->GetRandom();
+ }
+
+ // --- Uniform distribution: delta_t is the inverse event rate
+ else if (fTimeDist == TimeDist::Uniform) {
+ return 1.e9 / fEventRate; // rate is in 1/s, delta_t in ns
+ }
+
+ // --- Just in case: catch unknown distribution models
+ else
+ return 0.;
}
// ---------------------------------------------------------------------------
diff --git a/sim/response/base/CbmMCInputSet.h b/sim/response/base/CbmMCInputSet.h
index ec4960c39c..be725502d5 100644
--- a/sim/response/base/CbmMCInputSet.h
+++ b/sim/response/base/CbmMCInputSet.h
@@ -20,6 +20,8 @@
#include <set>
#include <utility>
+#include "Defs.h"
+
/** @class CbmMCInputSet
** @author Volker Friese <v.friese@gsi.de>
@@ -47,7 +49,7 @@ public:
/** @brief Constructor
** @param rate Event rate [1/s]. Must be positive.
**/
- CbmMCInputSet(Double_t rate);
+ CbmMCInputSet(cbm::sim::TimeDist dist, Double_t rate);
/** @brief Destructor **/
@@ -133,7 +135,7 @@ public:
/** @brief Event rate
** @value Event rate [1/s]
**/
- Double_t GetRate() const { return fRate; }
+ Double_t GetEventRate() const { return fEventRate; }
/** @brief register all input chains to the FairRootManager **/
@@ -141,11 +143,12 @@ public:
private:
- Double_t fRate; // Event rate [1/s]
- std::map<UInt_t, CbmMCInput*> fInputs; // Key is input ID
- std::map<UInt_t, CbmMCInput*>::iterator fInputHandle; // Handle for inputs
- std::set<TString> fBranches; // List of branch names
- TF1* fDeltaDist; // Probability distribution for delta(t)
+ cbm::sim::TimeDist fTimeDist = cbm::sim::TimeDist::Poisson; // Event time distribution model
+ Double_t fEventRate; // Event rate [1/s]
+ std::map<UInt_t, CbmMCInput*> fInputs; // Key is input ID
+ std::map<UInt_t, CbmMCInput*>::iterator fInputHandle; // Handle for inputs
+ std::set<TString> fBranches; // List of branch names
+ TF1* fDeltaDist; // Probability distribution for delta(t)
/** @brief Compare an input branch list with the reference branch list
diff --git a/sim/response/base/Defs.h b/sim/response/base/Defs.h
new file mode 100644
index 0000000000..0030057a5b
--- /dev/null
+++ b/sim/response/base/Defs.h
@@ -0,0 +1,37 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Volker Friese [committer] */
+
+/** @file CbmMCInput.h
+ ** @author Volker Friese <v.friese@gsi.de>
+ ** @date 20.06.2023
+ **
+ ** Definitions in namespace cbm::sim
+ **/
+
+#ifndef SIM_RESPONSE_BASE_DEFS_H
+#define SIM_RESPONSE_BASE_DEFS_H 1
+
+
+namespace cbm::sim
+{
+
+ /** @enum Simulation mode **/
+ enum class Mode
+ {
+ Timebased,
+ EventByEvent
+ };
+
+ /** @enum Time distribution model **/
+ enum class TimeDist
+ {
+ Poisson,
+ Uniform
+ };
+
+
+} // namespace cbm::sim
+
+
+#endif /* SIM_RESPONSE_BASE_DEFS_H */
--
GitLab