diff --git a/core/detectors/rich/CbmRichDigiMapManager.cxx b/core/detectors/rich/CbmRichDigiMapManager.cxx
index 087740cbdf0408f3674e2f401488262011db1dc9..5e264b0e5cd36f455694cf7df4d22de8c3f96067 100644
--- a/core/detectors/rich/CbmRichDigiMapManager.cxx
+++ b/core/detectors/rich/CbmRichDigiMapManager.cxx
@@ -88,10 +88,10 @@ void CbmRichDigiMapManager::Init()
case 0: {
// RICH detector
//std::cout<<"RICH"<<std::endl;
- const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
- const Double_t* curNodeTr = curMatrix->GetTranslation();
+ const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
+ const Double_t* curNodeTr = curMatrix->GetTranslation();
const TGeoMatrix* localMatrix = curNode->GetMatrix(); // local transformation of pixel in MAPMT volume
- string path = string(nodePath.Data());
+ string path = string(nodePath.Data());
size_t pmtInd = path.find_last_of("/");
if (string::npos == pmtInd) continue;
@@ -143,10 +143,10 @@ void CbmRichDigiMapManager::Init()
case 1: {
//std::cout<<"mRICH"<<std::endl;
- const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
- const Double_t* curNodeTr = curMatrix->GetTranslation();
+ const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
+ const Double_t* curNodeTr = curMatrix->GetTranslation();
const TGeoMatrix* localMatrix = curNode->GetMatrix(); // local transformation of pixel in MAPMT volume
- string path = string(nodePath.Data());
+ string path = string(nodePath.Data());
size_t pmtInd = path.find_last_of("/");
if (string::npos == pmtInd) continue;
@@ -322,9 +322,10 @@ CbmRichPmtData* CbmRichDigiMapManager::GetPmtDataById(Int_t id)
return it->second;
}
-vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t N, Bool_t direct, Bool_t diagonal)
+vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t n, Bool_t direct, Bool_t diagonal)
{
vector<Int_t> neighbourPixels;
+ if (n == 0) return neighbourPixels;
CbmRichPixelData* addressPixelData = GetPixelDataByAddress(address);
Int_t indX = addressPixelData->fPixelId % 8;
Int_t indY = addressPixelData->fPixelId / 8;
@@ -334,18 +335,18 @@ vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t N,
CbmRichPixelData* iPixelData = GetPixelDataByAddress(iAddr);
Int_t iIndX = iPixelData->fPixelId % 8;
Int_t iIndY = iPixelData->fPixelId / 8;
- if (direct && !diagonal && N == 1) { // direct neighbours
+ if (direct && !diagonal && n == 1) { // direct neighbours
if ((abs(iIndX - indX) + abs(iIndY - indY)) == 1) neighbourPixels.push_back(iAddr);
}
- else if (!direct && diagonal && N == 1) { // diagonal neighbours
+ else if (!direct && diagonal && n == 1) { // diagonal neighbours
if ((abs(iIndX - indX) == 1) && (abs(iIndY - indY) == 1)) neighbourPixels.push_back(iAddr);
}
else if (direct && diagonal) { // all neighbours in (2N+1)*(2N+1) grid
- if ((abs(iIndX - indX) <= N) && (abs(iIndY - indY) <= N)) neighbourPixels.push_back(iAddr);
+ if ((abs(iIndX - indX) <= n) && (abs(iIndY - indY) <= n)) neighbourPixels.push_back(iAddr);
}
else {
LOG(error) << "ERROR: Unrecogniced option in CbmRichDigiMapManager::GetNeighbourPixels " << endl
- << " N = " << N << " direct = " << direct << " diagonal = " << diagonal;
+ << " n = " << n << " direct = " << direct << " diagonal = " << diagonal;
}
}
return neighbourPixels;
diff --git a/core/detectors/rich/CbmRichDigiMapManager.h b/core/detectors/rich/CbmRichDigiMapManager.h
index 54a6323b279d6cbc37db3f1708e9acb57933d7e5..d4a9604081b12fb4b047ee46558af8ca4775cf99 100644
--- a/core/detectors/rich/CbmRichDigiMapManager.h
+++ b/core/detectors/rich/CbmRichDigiMapManager.h
@@ -68,28 +68,36 @@ public:
/*
* \brief Return the addresses of the neighbour pixels.
+ * \param address Pixel address
+ * \param n Size of the grid (2n+1)*(2n+1)
+ * \param direct return direct neighbours
+ * \param diagonal return diagonal neighbours
*/
std::vector<Int_t> GetNeighbourPixels(Int_t address, Int_t N, Bool_t direct = true, Bool_t diagonal = true);
/*
* \brief Return the addresses of the direct neighbour pixels.
+ * \param address Pixel address
*/
std::vector<Int_t> GetDirectNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, true, false); }
/*
* \brief Return the addresses of the diagonal neighbour pixels.
+ * \param address Pixel address
*/
std::vector<Int_t> GetDiagonalNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, false, true); }
/*
- * \brief Return the addresses of pixels in a (2N+1)*(2N+1) grid,
+ * \brief Return the addresses of pixels in a (2n+1)*(2n+1) grid,
* with the address pixel in the center of the grid.
* Addresses are limited to the same MAPMT as the input address.
* Needed for noise digis caused by charged particles.
+ * \param address Pixel address
+ * \param n Size of the grid (2n+1)*(2n+1)
*/
- std::vector<Int_t> GetNxNNeighbourPixels(Int_t address, Int_t N)
+ std::vector<Int_t> GetNxNNeighbourPixels(Int_t address, Int_t n)
{
- return GetNeighbourPixels(address, N, true, true);
+ return GetNeighbourPixels(address, n, true, true);
}
public:
diff --git a/macro/rich/geotest/run_digi.C b/macro/rich/geotest/run_digi.C
index 1ebc322414f6f43642a231961e2645d71dc0f969..85dacd37f4a928878c7531c00ead6769bd73c25f 100644
--- a/macro/rich/geotest/run_digi.C
+++ b/macro/rich/geotest/run_digi.C
@@ -27,8 +27,7 @@ void run_digi(const string& traFile, const string& parFile, const string& digiFi
run.SetMonitorFile("");
CbmRichDigitizer* richDigitizer = new CbmRichDigitizer();
- richDigitizer->SetMaxNofHitsPerPmtCut(65);
- richDigitizer->SetNoiseDigiRate(0.);
+ richDigitizer->SetEventNoiseRate(0.);
run.SetDigitizer(ECbmModuleId::kRich, richDigitizer, "RichPoint", true);
run.Run(nEvents);
diff --git a/macro/rich/run/run_digi.C b/macro/rich/run/run_digi.C
index 8b91747631f48b9af7d912882d59571c5c3509f4..20093cb2628cdba4385ce46dbf356bbdec73eeed 100644
--- a/macro/rich/run/run_digi.C
+++ b/macro/rich/run/run_digi.C
@@ -31,7 +31,6 @@ void run_digi(const string& traFile = "/Users/slebedev/Development/cbm/data/sim
run.SetMonitorFile("");
CbmRichDigitizer* richDigitizer = new CbmRichDigitizer();
- richDigitizer->SetMaxNofHitsPerPmtCut(65);
run.SetDigitizer(ECbmModuleId::kRich, richDigitizer, "RichPoint", true);
run.Run(nEvents);
diff --git a/reco/detectors/rich/CMakeLists.txt b/reco/detectors/rich/CMakeLists.txt
index 7fecbf83e010114b46eaceebb273453eff64791c..becc3e386c136430fc70b1efdd8ecc87c5abab77 100644
--- a/reco/detectors/rich/CMakeLists.txt
+++ b/reco/detectors/rich/CMakeLists.txt
@@ -34,6 +34,7 @@ set(SRCS
qa/CbmRichGeoTestOpt.cxx
qa/CbmRichRecoQa.cxx
qa/CbmL1RichRingQa.cxx
+ qa/CbmRichDigiQa.cxx
unpack/CbmRichUnpackAlgo.cxx
unpack/CbmRichUnpackAlgo2022.cxx
diff --git a/reco/detectors/rich/CbmRichRecoLinkDef.h b/reco/detectors/rich/CbmRichRecoLinkDef.h
index 8739d2c8657490ca67fcc988c48c37ab0dd8b886..a832c0dd25658ce645712bbbd8b6586ee53ffd7c 100644
--- a/reco/detectors/rich/CbmRichRecoLinkDef.h
+++ b/reco/detectors/rich/CbmRichRecoLinkDef.h
@@ -32,6 +32,7 @@
#pragma link C++ class CbmRichRecoQa + ;
#pragma link C++ class CbmRichRecoTbQa + ;
#pragma link C++ class CbmL1RichRingQa + ;
+#pragma link C++ class CbmRichDigiQa + ;
//unpack
#pragma link C++ class CbmRichUnpackAlgo + ;
diff --git a/reco/detectors/rich/qa/CbmRichDigiQa.cxx b/reco/detectors/rich/qa/CbmRichDigiQa.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4524ad7b3695664afa90a1f0ca9adcf15e85b9c4
--- /dev/null
+++ b/reco/detectors/rich/qa/CbmRichDigiQa.cxx
@@ -0,0 +1,123 @@
+/* Copyright (C) 2023-2023 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+#include "CbmRichDigiQa.h"
+
+#include "CbmDigiManager.h"
+#include "CbmHistManager.h"
+#include "CbmRichDetectorData.h"
+#include "CbmRichDigi.h"
+#include "CbmRichDigiMapManager.h"
+
+#include <FairRootManager.h>
+#include <Logger.h>
+
+#include <TDirectory.h>
+#include <TFile.h>
+
+InitStatus CbmRichDigiQa::Init()
+{
+ CbmRichDigiMapManager::GetInstance();
+
+ fDigiMan = CbmDigiManager::Instance();
+ fDigiMan->Init();
+ if (!fDigiMan->IsPresent(ECbmModuleId::kRich)) LOG(fatal) << "CbmRichDigiQa::Init: No RichDigi array!";
+
+ InitHistograms();
+
+ return kSUCCESS;
+}
+
+void CbmRichDigiQa::InitHistograms()
+{
+ fHM = new CbmHistManager();
+
+ fHM->Create1<TH1D>("fhDigiDt", "same address consecutive digis time diff;dt [ns];entries", 1001, -0.5, 10000.5);
+
+ fHM->Create1<TH1D>("fhDigiDtEdge", "same address consecutive digis time diff;dt [ns];entries", 101, -0.5, 100.5);
+
+ fHM->Create1<TH1D>("fhDigiDistH", "horizontal distance of digis same pmt;DIndX [in pixels];entries", 11, -0.5, 10.5);
+
+ fHM->Create1<TH1D>("fhDigiDistV", "vertical distance of digis same pmt;DIndY [in pixels];entries", 11, -0.5, 10.5);
+
+ fHM->Create2<TH2D>("fhDigiNeighbours", "digi neighbours same pmt;DIndX [in pixels];DIndY [in pixels];entries", 5,
+ -2.5, 2.5, 5, -2.5, 2.5);
+}
+
+void CbmRichDigiQa::Exec(Option_t* /*option*/)
+{
+ LOG(debug) << GetName() << " Event " << fEventNum;
+ fEventNum++;
+
+ Int_t nofDigis = fDigiMan->GetNofDigis(ECbmModuleId::kRich);
+ for (Int_t i = 0; i < nofDigis; i++) {
+ const CbmRichDigi* digi = fDigiMan->Get<CbmRichDigi>(i);
+ if (!digi) continue;
+ Int_t pmtId = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(digi->GetAddress())->fPmtId;
+ if (fToTLimitLow > 0. && fToTLimitHigh > 0.) {
+ if (digi->GetToT() < fToTLimitHigh && digi->GetToT() > fToTLimitLow) {
+ fFiredTimes[digi->GetAddress()].push_back(digi->GetTime());
+ fPmtDigisTimeAddress[pmtId].push_back(std::make_pair(digi->GetTime(), digi->GetAddress()));
+ }
+ }
+ else {
+ fFiredTimes[digi->GetAddress()].push_back(digi->GetTime());
+ fPmtDigisTimeAddress[pmtId].push_back(std::make_pair(digi->GetTime(), digi->GetAddress()));
+ }
+ }
+
+ // Fill deadtimes
+ for (auto& [address, times] : fFiredTimes) {
+ std::sort(times.begin(), times.end());
+ for (auto it = times.begin(); it != times.end(); ++it) {
+ if (*it < 0.)
+ LOG(warn) << "Digi time < 0.: "
+ << "address: " << address << " time: " << *it << " ns";
+
+ if (it != times.begin()) {
+ Double_t dt = *it - *std::prev(it);
+ fHM->H1("fhDigiDt")->Fill(dt);
+ if (dt < 100.5) { fHM->H1("fhDigiDtEdge")->Fill(dt); }
+ }
+ }
+ }
+
+ // Fill crosstalk digis
+ for (auto& [pmt, digis] : fPmtDigisTimeAddress) {
+ std::sort(digis.begin(), digis.end());
+ for (auto it0 = digis.begin(); it0 != digis.end(); ++it0) {
+ Int_t indX0 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it0->second)->fPixelId % 8;
+ Int_t indY0 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it0->second)->fPixelId / 8;
+ for (auto it1 = it0 + 1; it1 != digis.end(); ++it1) {
+ if (abs(it0->first - it1->first) < fCrosstalkTimeLimit) {
+ Int_t indX1 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it1->second)->fPixelId % 8;
+ Int_t indY1 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it1->second)->fPixelId / 8;
+ if (indX0 == indX1 && indY0 == indY1) continue;
+ if (indY0 == indY1) fHM->H1("fhDigiDistH")->Fill(indX1 - indX0);
+ if (indX0 == indX1) fHM->H1("fhDigiDistV")->Fill(indY1 - indY0);
+ if (abs(indX0 - indX1) > 1 || abs(indY0 - indY1) > 1) continue;
+ fHM->H2("fhDigiNeighbours")->Fill(indX1 - indX0, indY0 - indY1);
+ }
+ else if (it1->first - it0->first > fCrosstalkTimeLimit) {
+ break;
+ }
+ }
+ }
+ }
+
+ fFiredTimes.clear();
+ fPmtDigisTimeAddress.clear();
+}
+
+void CbmRichDigiQa::Finish()
+{
+ TDirectory* oldir = gDirectory;
+ TFile* outFile = FairRootManager::Instance()->GetOutFile();
+ if (outFile != nullptr) {
+ outFile->mkdir(GetName());
+ outFile->cd(GetName());
+ fHM->WriteToFile();
+ }
+ gDirectory->cd(oldir->GetPath());
+}
diff --git a/reco/detectors/rich/qa/CbmRichDigiQa.h b/reco/detectors/rich/qa/CbmRichDigiQa.h
new file mode 100644
index 0000000000000000000000000000000000000000..25822a5c28a7266e4a63c1fbf2da77b5be8541f9
--- /dev/null
+++ b/reco/detectors/rich/qa/CbmRichDigiQa.h
@@ -0,0 +1,94 @@
+/* Copyright (C) 2023-2023 UGiessen, Giessen
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Martin Beyer [committer] */
+
+/**
+* \file CbmRichDigiQa.h
+* \author M.Beyer
+* \date 2023
+**/
+
+#ifndef CBMRICHDIGIQA_H
+#define CBMRICHDIGIQA_H
+
+#include "FairTask.h"
+
+#include <map>
+#include <vector>
+
+class CbmDigiManager;
+class CbmRichDigi;
+class CbmHistManager;
+class TClonesArray;
+
+/**
+* \class CbmRichDigiQa
+*
+* \brief Class for pixel deadtime and crosstalk Qa.
+*
+* \author M.Beyer
+* \date 2023
+**/
+class CbmRichDigiQa : public FairTask {
+public:
+ /** Default constructor */
+ CbmRichDigiQa() : FairTask("CbmRichDigiQa") {};
+
+ /** Destructor */
+ ~CbmRichDigiQa() = default;
+
+ /** Copy constructor (disabled) */
+ CbmRichDigiQa(const CbmRichDigiQa&) = delete;
+
+ /** Assignment operator (disabled) */
+ CbmRichDigiQa operator=(const CbmRichDigiQa&) = delete;
+
+ /** Inherited from FairTask */
+ InitStatus Init();
+
+ /** Inherited from FairTask */
+ void Exec(Option_t* option);
+
+ /** Inherited from FairTask */
+ void Finish();
+
+ /** Initialize histogram manager */
+ void InitHistograms();
+
+ /**
+ * \brief Set ToT cut
+ * \param[in] low lower limit
+ * \param[in] high upper limit
+ */
+ void SetToTLimits(Double_t low, Double_t high)
+ {
+ fToTLimitLow = low;
+ fToTLimitHigh = high;
+ }
+
+ /**
+ * \brief Set time limit in which a neighbour digi is considered crosstalk
+ * \param[in] limit absolue time limit
+ */
+ void SetCrosstalkTimeLimit(Double_t limit) { fCrosstalkTimeLimit = limit; }
+
+private:
+ int fEventNum {};
+
+ CbmDigiManager* fDigiMan {nullptr};
+ CbmHistManager* fHM {nullptr};
+
+ TClonesArray* fRichPoints {nullptr};
+
+ std::map<Int_t, std::vector<Double_t>> fFiredTimes {}; // Store times of fired digis for each pixel in Ts
+ std::map<Int_t, std::vector<std::pair<Double_t, Int_t>>>
+ fPmtDigisTimeAddress {}; // Store times and addresses of fired digis for each pmt in Ts
+
+ Double_t fCrosstalkTimeLimit {5.}; // Time in which a neighbour digi is considered crosstalk
+ Double_t fToTLimitLow {-1.}; // Lower limit for ToT cut
+ Double_t fToTLimitHigh {-1.}; // Upper limit for ToT cut
+
+ ClassDef(CbmRichDigiQa, 1)
+};
+
+#endif
diff --git a/sim/detectors/rich/CbmRichDigitizer.cxx b/sim/detectors/rich/CbmRichDigitizer.cxx
index 6b6761cd6ed0154b02a3b7d14db50bfcc066eefc..dedeebf2c3e0175cf588d975c4994417ff097ecb 100644
--- a/sim/detectors/rich/CbmRichDigitizer.cxx
+++ b/sim/detectors/rich/CbmRichDigitizer.cxx
@@ -1,6 +1,6 @@
-/* Copyright (C) 2015-2020 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
SPDX-License-Identifier: GPL-3.0-only
- Authors: Andrey Lebedev [committer], Volker Friese, Florian Uhlig, Semen Lebedev */
+ Authors: Martin Beyer, Andrey Lebedev [committer], Volker Friese, Florian Uhlig, Semen Lebedev */
/**
* \file CbmRichDigitizer.cxx
@@ -14,64 +14,39 @@
#include "CbmLink.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
-#include "CbmRichDetectorData.h" // for CbmRichPmtData, CbmRichPixelData
+#include "CbmRichDetectorData.h"
#include "CbmRichDigi.h"
#include "CbmRichDigiMapManager.h"
#include "CbmRichGeoManager.h"
#include "CbmRichPoint.h"
-#include "FairEventHeader.h"
-#include "FairRunAna.h"
-#include "FairRunSim.h"
+#include <FairEventHeader.h>
+#include <FairRunAna.h>
+#include <FairRunSim.h>
#include <Logger.h>
-#include "TClonesArray.h"
-#include "TGeoManager.h"
-#include "TRandom.h"
-#include "TStopwatch.h"
+#include <TClonesArray.h>
+#include <TGeoManager.h>
+#include <TMath.h>
+#include <TRandom.h>
+#include <TStopwatch.h>
#include <cassert>
#include <iomanip>
#include <iostream>
-using namespace std;
-
-CbmRichDigitizer::CbmRichDigitizer()
- : CbmDigitize<CbmRichDigi>("RichDigitize")
- , fEventNum(0)
- , fRichPoints(NULL)
- , fRichDigis(NULL)
- , fRichDigiMatches(nullptr)
- , fMcTracks(NULL)
- , fNofPoints(0.)
- , fNofDigis(0.)
- , fTimeTot(0.)
- , fPmt()
- , fCrossTalkProbability(0.02)
- , fNoiseDigiRate(5.)
- , fDetectorType(CbmRichPmtTypeCosy17NoWls)
- , fMaxNofHitsPerPmtCut(65)
- , fDataMap()
- , fTimeResolution(1.)
- , fDarkRatePerPixel(1000)
- , fPixelDeadTime(30.)
- , fFiredPixelsMap()
- , fDoZShift(true)
-{
-}
-
-CbmRichDigitizer::~CbmRichDigitizer() {}
-
+using std::fixed;
+using std::left;
+using std::right;
+using std::setprecision;
+using std::setw;
InitStatus CbmRichDigitizer::Init()
{
-
- // Screen output
std::cout << std::endl;
LOG(info) << "==========================================================";
LOG(info) << GetName() << ": Initialisation";
- //FairTask* daq = FairRun::Instance()->GetTask("Daq");
if (!fEventMode) { LOG(info) << "CbmRichDigitizer uses TimeBased mode."; }
else {
LOG(info) << "CbmRichDigitizer uses Events mode.";
@@ -84,11 +59,11 @@ InitStatus CbmRichDigitizer::Init()
CbmRichDigiMapManager::GetInstance();
CbmRichGeoManager::GetInstance();
- fRichPoints = (TClonesArray*) manager->GetObject("RichPoint");
- if (NULL == fRichPoints) { Fatal("CbmRichDigitizer::Init", "No RichPoint array!"); }
+ fRichPoints = static_cast<TClonesArray*>(manager->GetObject("RichPoint"));
+ if (!fRichPoints) LOG(fatal) << "CbmRichDigitizer::Init: No RichPoint array!";
- fMcTracks = (TClonesArray*) manager->GetObject("MCTrack");
- if (NULL == fMcTracks) { Fatal("CbmRichDigitizer::Init", "No MCTrack array!"); }
+ fMcTracks = static_cast<TClonesArray*>(manager->GetObject("MCTrack"));
+ if (!fMcTracks) LOG(fatal) << "CbmRichDigitizer::Init: No MCTrack array!";
RegisterOutput();
@@ -105,8 +80,8 @@ InitStatus CbmRichDigitizer::Init()
}
LOG(info) << GetName() << ": Initialisation successful";
-
- fFiredPixelsMap.clear();
+ LOG(info) << "==========================================================";
+ std::cout << std::endl;
return kSUCCESS;
}
@@ -116,31 +91,21 @@ void CbmRichDigitizer::Exec(Option_t* /*option*/)
TStopwatch timer;
timer.Start();
- fEventNum++;
- LOG(debug) << fName << ": Event " << fEventNum;
-
-
- // Clear data map
- for (auto& data : fDataMap) {
- if (data.second.first) delete data.second.first; // digi
- if (data.second.second) delete data.second.second; // match
- }
- fDataMap.clear();
-
Double_t oldEventTime = fCurrentEventTime;
GetEventInfo();
- LOG(debug) << fName << ": EventNumber: " << fCurrentEvent << ", EventTime: " << fCurrentEventTime;
-
- if (fProduceNoise) GenerateNoiseBetweenEvents(oldEventTime, fCurrentEventTime);
-
Int_t nPoints = ProcessMcEvent();
- Int_t nDigis = AddDigisToOutputArray();
+ Double_t tNoiseStart = fEventNum ? oldEventTime : fRunStartTime;
+ if (fProduceNoise) AddDarkRateNoise(tNoiseStart, fCurrentEventTime);
+ if (fEventMode) fPixelsLastFiredTime.clear();
+ Double_t processTime = fEventMode ? -1. : fCurrentEventTime;
+ Int_t nDigis = ProcessBuffers(processTime);
// --- Statistics
timer.Stop();
+ fEventNum++;
fNofPoints += nPoints;
fNofDigis += nDigis;
fTimeTot += timer.RealTime();
@@ -151,9 +116,36 @@ void CbmRichDigitizer::Exec(Option_t* /*option*/)
<< ". Exec time " << setprecision(6) << timer.RealTime() << " s.";
}
+void CbmRichDigitizer::Finish()
+{
+ TStopwatch timer;
+ if (!fEventMode) {
+ std::cout << std::endl;
+ LOG(info) << GetName() << ": Finish run";
+ LOG(info) << GetName() << ": Processing analogue buffers";
+ timer.Start();
+ Int_t nDigis = ProcessBuffers(-1.);
+ timer.Stop();
+ fTimeTot += timer.RealTime();
+ LOG(info) << GetName() << ": " << nDigis << " digis created and sent to DAQ";
+ fNofDigis += nDigis;
+ }
+
+ std::cout << std::endl;
+ LOG(info) << "=====================================";
+ LOG(info) << GetName() << ": Run summary";
+ LOG(info) << "Events processed : " << fEventNum;
+ LOG(info) << "RichPoint / event : " << setprecision(1) << fNofPoints / Double_t(fEventNum);
+ LOG(info) << "RichDigi / event : " << fNofDigis / Double_t(fEventNum);
+ LOG(info) << "Digis per point : " << setprecision(6) << fNofDigis / Double_t(fNofPoints);
+ LOG(info) << "Real time per event : " << fTimeTot / Double_t(fEventNum) << " s";
+ LOG(info) << "=====================================";
+}
+
Int_t CbmRichDigitizer::ProcessMcEvent()
{
Int_t nofRichPoints = fRichPoints->GetEntriesFast();
+
LOG(debug) << fName << ": EventNum:" << fCurrentEvent << " InputNum:" << fCurrentInput
<< " EventTime:" << fCurrentEventTime << " nofRichPoints:" << nofRichPoints;
@@ -161,30 +153,12 @@ Int_t CbmRichDigitizer::ProcessMcEvent()
CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(j);
ProcessPoint(point, j, fCurrentMCEntry, fCurrentInput);
}
- // cout << "nofDigis:" << fRichDigis->GetEntriesFast() << endl;
- AddNoiseDigis(fCurrentMCEntry, fCurrentInput);
+ AddEventNoise(fCurrentMCEntry, fCurrentInput);
return nofRichPoints;
}
-void CbmRichDigitizer::GenerateNoiseBetweenEvents(Double_t oldEventTime, Double_t newEventTime)
-{
- Int_t nofPixels = CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size();
- Double_t dT = newEventTime - oldEventTime;
- Double_t nofNoisePixels = nofPixels * dT * (fDarkRatePerPixel / 1.e9);
-
- LOG(debug) << "CbmRichDigitizer::GenerateNoiseBetweenEvents deltaTime:" << dT << " nofPixels:" << nofPixels
- << " nofNoisePixels:" << nofNoisePixels;
-
- for (Int_t j = 0; j < nofNoisePixels; j++) {
- Int_t address = CbmRichDigiMapManager::GetInstance().GetRandomPixelAddress();
- CbmLink link(1., -1, -1, -1);
- Double_t time = gRandom->Uniform(oldEventTime, newEventTime);
- AddDigi(address, time, link);
- }
-}
-
void CbmRichDigitizer::ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t eventNum, Int_t inputNum)
{
// LOG(info) << "ProcessPoint " << pointId;
@@ -193,15 +167,16 @@ void CbmRichDigitizer::ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t ev
Double_t zNew = point->GetZ();
if (fDoZShift) zNew = zNew - 0.001;
gGeoManager->FindNode(point->GetX(), point->GetY(), zNew);
- string path(gGeoManager->GetPath());
+ std::string path(gGeoManager->GetPath());
Int_t address = CbmRichDigiMapManager::GetInstance().GetPixelAddressByPath(path);
-
+ if (address < 0) return;
Int_t trackId = point->GetTrackID();
//LOG(info) << "address:" << address << " path:" << path << " trackId:" << trackId << "X:" << point->GetX() << " Y:" << point->GetY() << " Z:" << zNew;
if (trackId < 0) return;
CbmMCTrack* p = (CbmMCTrack*) fMcTracks->At(trackId);
if (p == NULL) return;
Int_t gcode = TMath::Abs(p->GetPdgCode());
+ Double_t charge = p->GetCharge();
CbmLink link(1., pointId, eventNum, inputNum);
@@ -219,165 +194,143 @@ void CbmRichDigitizer::ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t ev
isDetected = true;
}
- Double_t time = fCurrentEventTime + point->GetTime();
- Double_t deltaT = gRandom->Gaus(0., fTimeResolution);
- time += deltaT;
- if (isDetected) { AddDigi(address, time, link); }
+ Double_t time = fCurrentEventTime + point->GetTime();
+ if (isDetected) {
+ AddSignalToBuffer(address, time, link);
+ if (gcode == 50000050 && address > 0) AddCrossTalk(address, time, link);
+ if (charge != 0. && address > 0) AddChargedParticleCluster(address, time, eventNum, inputNum);
+ }
}
-
-void CbmRichDigitizer::Finish()
+void CbmRichDigitizer::AddSignalToBuffer(Int_t address, Double_t time, const CbmLink& link)
{
- std::cout << std::endl;
- LOG(info) << "=====================================";
- LOG(info) << GetName() << ": Run summary";
- LOG(info) << "Events processed : " << fEventNum;
- LOG(info) << "RichPoint / event : " << setprecision(1) << fNofPoints / Double_t(fEventNum);
- LOG(info) << "RichDigi / event : " << fNofDigis / Double_t(fEventNum);
- LOG(info) << "Digis per point : " << setprecision(6) << fNofDigis / fNofPoints;
- LOG(info) << "Real time per event : " << fTimeTot / Double_t(fEventNum) << " s";
- LOG(info) << "=====================================";
+ fSignalBuffer[address].emplace_back(std::make_pair(time, new CbmLink(link)));
}
-void CbmRichDigitizer::AddCrossTalkDigis(CbmRichDigi* digi)
+void CbmRichDigitizer::AddCrossTalk(Int_t address, Double_t time, const CbmLink& link)
{
- Bool_t crosstalkDigiDetected = false;
- vector<Int_t> directPixels = CbmRichDigiMapManager::GetInstance().GetDirectNeighbourPixels(digi->GetAddress());
- vector<Int_t> diagonalPixels = CbmRichDigiMapManager::GetInstance().GetDiagonalNeighbourPixels(digi->GetAddress());
-
- for (UInt_t i = 0; i < directPixels.size(); i++) {
- if (gRandom->Rndm() < fCrossTalkProbability && !crosstalkDigiDetected) {
- //FindRichHitPositionMAPMT(0., x + r, y, xHit, yHit, pmtID);
- crosstalkDigiDetected = true;
- break;
- }
- }
-
- if (!crosstalkDigiDetected) {
- for (UInt_t i = 0; i < diagonalPixels.size(); i++) {
- if (gRandom->Rndm() < fCrossTalkProbability / 4. && !crosstalkDigiDetected) {
- //FindRichHitPositionMAPMT(0., x + r, y, xHit, yHit, pmtID);
- crosstalkDigiDetected = true;
- break;
- }
+ std::vector<Int_t> directPixels = CbmRichDigiMapManager::GetInstance().GetDirectNeighbourPixels(address);
+ std::vector<Int_t> diagonalPixels = CbmRichDigiMapManager::GetInstance().GetDiagonalNeighbourPixels(address);
+ /* clang-format off */
+ Double_t crosstalkCreatedProbability = 1 - pow(1 - fCrossTalkProbability, directPixels.size()) *
+ pow(1 - fCrossTalkProbability / 4., diagonalPixels.size());
+
+ Int_t crosstalkAddress = -1;
+ if (gRandom->Uniform() < crosstalkCreatedProbability) {
+ Double_t thresholdDirectDiagonal = directPixels.size() / float(directPixels.size() + diagonalPixels.size() / 4.);
+ if (gRandom->Uniform() < thresholdDirectDiagonal) {
+ crosstalkAddress = directPixels[gRandom->Integer(directPixels.size())];
+ } else {
+ crosstalkAddress = diagonalPixels[gRandom->Integer(diagonalPixels.size())];
}
}
-
- if (crosstalkDigiDetected) {
- //AddDigi(posHit, posHitErr, RichDetID, pmtID, ampl, pointInd);
- //fNofCrossTalkDigis++;
- }
+ /* clang-format on */
+ if (crosstalkAddress != -1) AddSignalToBuffer(crosstalkAddress, time, link);
}
-Int_t CbmRichDigitizer::AddDigisToOutputArray()
+void CbmRichDigitizer::AddChargedParticleCluster(Int_t address, Double_t time, Int_t eventNum, Int_t inputNum)
{
- Int_t nofDigis = 0;
-
- // fMaxNofHitsPerPmtCut is a maximum number of hits which can be registered per PMT per event.
- // If more then the whole PMT is skipped
- map<Int_t, Int_t> digisPerPmtMap;
- for (auto const& dm : fDataMap) {
- CbmRichPixelData* pixelData =
- CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(dm.second.first->GetAddress());
- if (nullptr == pixelData) continue;
- Int_t pmtId = pixelData->fPmtId;
- digisPerPmtMap[pmtId]++;
- }
-
- Int_t nofSkippedPmts = 0;
- for (auto const& dm : digisPerPmtMap) {
- if (dm.second > fMaxNofHitsPerPmtCut) nofSkippedPmts++;
- }
-
- Int_t nofSkippedDigis = 0;
- for (auto const& dm : fDataMap) {
- CbmRichPixelData* pixelData =
- CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(dm.second.first->GetAddress());
- if (nullptr == pixelData) continue;
- if (digisPerPmtMap[pixelData->fPmtId] > fMaxNofHitsPerPmtCut) {
- nofSkippedDigis++;
- continue;
+ Int_t sourcePixelId = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(address)->fPixelId;
+ std::vector<Int_t> neighbourAddresses =
+ CbmRichDigiMapManager::GetInstance().GetNxNNeighbourPixels(address, fClusterSize);
+ for (UInt_t i = 0; i < neighbourAddresses.size(); i++) {
+ Int_t iPixelId = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(neighbourAddresses[i])->fPixelId;
+ Double_t d = TMath::Sqrt((sourcePixelId % 8 - iPixelId % 8) * (sourcePixelId % 8 - iPixelId % 8)
+ + (sourcePixelId / 8 - iPixelId / 8) * (sourcePixelId / 8 - iPixelId / 8));
+ if (gRandom->Uniform(1.) < fClusterSignalProbability / d) {
+ CbmLink link(1., -1, eventNum, inputNum);
+ AddSignalToBuffer(neighbourAddresses[i], time, link);
}
-
- CbmRichDigi* digi = new CbmRichDigi();
- digi->SetAddress(dm.second.first->GetAddress());
- CbmMatch* digiMatch = new CbmMatch(*dm.second.second);
- digi->SetTime(dm.second.first->GetTime());
- //SendDigi(digi, digiMatch);
- if (fCreateMatches) SendData(digi->GetTime(), digi, digiMatch);
- else
- SendData(digi->GetTime(), digi);
- nofDigis++;
- } //# digis in map
-
- LOG(debug) << "Number of RICH digis before skip: " << fDataMap.size();
- LOG(debug) << "Nof skipped PMTs:" << nofSkippedPmts;
- LOG(debug) << "Nof skipped digis:" << nofSkippedDigis;
- LOG(debug) << "Number of RICH digis: " << nofDigis;
- return nofDigis;
+ }
}
-void CbmRichDigitizer::AddNoiseDigis(Int_t eventNum, Int_t inputNum)
+void CbmRichDigitizer::AddEventNoise(Int_t eventNum, Int_t inputNum)
{
Int_t nofPixels = CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size();
- Double_t dT = 50.; //ns
Double_t nofRichPoints = fRichPoints->GetEntriesFast();
- Double_t nofNoiseDigis = nofRichPoints * nofPixels * dT * (fNoiseDigiRate / 1.e9);
+ Double_t nofNoiseDigis = nofRichPoints * nofPixels * fEventNoiseInterval * (fEventNoiseRate / 1.e9);
+
+ LOG(debug) << "CbmRichDigitizer::AddEventNoise NofAllPixels:" << nofPixels << " nofNoiseDigis:" << nofNoiseDigis;
- LOG(debug) << "CbmRichDigitizer::AddNoiseDigis NofAllPixels:" << nofPixels << " nofNoiseDigis:" << nofNoiseDigis;
for (Int_t j = 0; j < nofNoiseDigis; j++) {
Int_t address = CbmRichDigiMapManager::GetInstance().GetRandomPixelAddress();
CbmLink link(1., -1, eventNum, inputNum);
- Double_t time = fCurrentEventTime + gRandom->Uniform(0, dT);
- AddDigi(address, time, link);
+ Double_t time = fCurrentEventTime + gRandom->Uniform(0, fEventNoiseInterval);
+ AddSignalToBuffer(address, time, link);
}
}
-void CbmRichDigitizer::AddDigi(Int_t address, Double_t time, const CbmLink& link)
+void CbmRichDigitizer::AddDarkRateNoise(Double_t oldEventTime, Double_t newEventTime)
{
- Bool_t wasFired = fFiredPixelsMap.count(address) > 0;
- Bool_t isDetected = true;
- if (!fEventMode && wasFired) {
- Double_t lastFiredTime = fFiredPixelsMap[address];
- Double_t dt = std::fabs(time - lastFiredTime);
- if (dt < fPixelDeadTime) {
- isDetected = false;
- // LOG(info) << "CbmRichDigitizer::AddDigi pixel NOT registered: address:" << address << " cur-last=dT: "<< time << "-"
- // << lastFiredTime << "=" << dt << " fPixelDeadTime:" << fPixelDeadTime
- // << " fFiredPixelsMap.size():" << fFiredPixelsMap.size() << " link: " << link.GetIndex();
- }
- else {
- isDetected = true;
- // LOG(info) << "CbmRichDigitizer::AddDigi pixel registered: address:" << address << " cur-last=dT: "<< time << "-"
- // << lastFiredTime << "=" << dt << " fPixelDeadTime:" << fPixelDeadTime
- // << " fFiredPixelsMap.size():" << fFiredPixelsMap.size() << " link: " << link.GetIndex();
- }
- }
- else {
- isDetected = true;
- // LOG(info) << "CbmRichDigitizer::AddDigi pixel was not fired before: address:" << address << " time:" << time
- // << " fFiredPixelsMap.size():" << fFiredPixelsMap.size() << " link: " << link.GetIndex();
+ Int_t nofPixels = CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size();
+ Double_t dT = newEventTime - oldEventTime;
+ Double_t nofNoisePixels = nofPixels * dT * (fDarkRatePerPixel / 1.e9);
+
+ LOG(debug) << "CbmRichDigitizer::AddDarkRateNoise deltaTime:" << dT << " nofPixels:" << nofPixels
+ << " nofNoisePixels:" << nofNoisePixels;
+
+ for (Int_t j = 0; j < nofNoisePixels; j++) {
+ Int_t address = CbmRichDigiMapManager::GetInstance().GetRandomPixelAddress();
+ CbmLink link(1., -1, -1, -1);
+ Double_t time = gRandom->Uniform(oldEventTime, newEventTime);
+ AddSignalToBuffer(address, time, link);
}
+}
- if (isDetected) {
- auto it = fDataMap.find(address);
- if (it == fDataMap.end()) {
- CbmRichDigi* digi = new CbmRichDigi();
- digi->SetAddress(address);
- digi->SetTime(time);
- CbmMatch* match = new CbmMatch();
- match->AddLink(link);
- pair<CbmRichDigi*, CbmMatch*> value(digi, match);
- fDataMap.insert(pair<Int_t, pair<CbmRichDigi*, CbmMatch*>>(address, value));
+Int_t CbmRichDigitizer::ProcessBuffers(Double_t processTime)
+{
+ Double_t maxNewDigiTime = processTime - fPixelDeadTime;
+ Int_t nofDigis {};
+
+ for (auto& dm : fSignalBuffer) {
+ std::sort(
+ dm.second.begin(), dm.second.end(),
+ [](const std::pair<Double_t, CbmLink*>& a, const std::pair<Double_t, CbmLink*>& b) { return a.first < b.first; });
+
+ CbmRichDigi* digi = nullptr;
+ CbmMatch* digiMatch = nullptr;
+ auto it = dm.second.begin();
+ for (; it != dm.second.end(); ++it) {
+ Double_t signalTime = it->first;
+ if (processTime >= 0. && signalTime > processTime) break;
+ CbmLink* link = it->second;
+ Bool_t createsDigi = true;
+
+ auto itFpm = fPixelsLastFiredTime.find(dm.first);
+ if (itFpm != fPixelsLastFiredTime.end()) {
+ Double_t lastFiredTime = itFpm->second;
+ if (signalTime - lastFiredTime < fPixelDeadTime) createsDigi = false;
+ }
+
+ if (createsDigi) {
+ if (digi) { // Add previous digi if exists
+ fCreateMatches ? SendData(digi->GetTime(), digi, digiMatch) : SendData(digi->GetTime(), digi);
+ nofDigis++;
+ digi = nullptr; // Needed because of potential break statement
+ digiMatch = nullptr;
+ }
+ Double_t digiTime = signalTime + gRandom->Gaus(0., fTimeResolution);
+ if (processTime >= 0. && digiTime > maxNewDigiTime) break;
+ digi = new CbmRichDigi();
+ digi->SetAddress(dm.first);
+ digi->SetTime(digiTime);
+ digiMatch = new CbmMatch();
+ digiMatch->AddLink(*link);
+ fPixelsLastFiredTime[dm.first] = signalTime; // Deadtime is from signal to signal
+ }
+ else {
+ if (digi && digiMatch) digiMatch->AddLink(*link);
+ }
+ delete it->second;
}
- else {
- CbmMatch* match = it->second.second;
- match->AddLink(link);
+ dm.second.erase(dm.second.begin(), it);
+ if (digi) { // Add last digi if exists
+ fCreateMatches ? SendData(digi->GetTime(), digi, digiMatch) : SendData(digi->GetTime(), digi);
+ nofDigis++;
}
}
- fFiredPixelsMap[address] = time;
+ return nofDigis;
}
-
ClassImp(CbmRichDigitizer)
diff --git a/sim/detectors/rich/CbmRichDigitizer.h b/sim/detectors/rich/CbmRichDigitizer.h
index 9c0bf86b7fa7a0d60a52b7cfed983b54b90152e5..6283766327bbe9157ecaa1ca91cf6d9534cad87b 100644
--- a/sim/detectors/rich/CbmRichDigitizer.h
+++ b/sim/detectors/rich/CbmRichDigitizer.h
@@ -1,6 +1,6 @@
-/* Copyright (C) 2015-2020 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
SPDX-License-Identifier: GPL-3.0-only
- Authors: Semen Lebedev, Andrey Lebedev [committer] */
+ Authors: Semen Lebedev, Andrey Lebedev [committer], Martin Beyer */
/**
* \file CbmRichDigitizer.h
@@ -23,12 +23,9 @@
#include <map>
class TClonesArray;
-//class CbmRichDigi;
class CbmRichPoint;
class CbmLink;
-using namespace std;
-
/**
* \class CbmRichDigitizer
@@ -40,79 +37,82 @@ using namespace std;
**/
class CbmRichDigitizer : public CbmDigitize<CbmRichDigi> {
public:
- /**
- * \brief Default constructor.
- */
- CbmRichDigitizer();
-
- /**
- * \brief Destructor.
- */
- virtual ~CbmRichDigitizer();
+ /** Default constructor */
+ CbmRichDigitizer() : CbmDigitize<CbmRichDigi>("RichDigitize") {};
+ /** Destructor */
+ virtual ~CbmRichDigitizer() = default;
- /** @brief Detector system ID
- ** @return kRich
- **/
- ECbmModuleId GetSystemId() const { return ECbmModuleId::kRich; }
+ /** Copy constructor (disabled) */
+ CbmRichDigitizer(const CbmRichDigitizer&) = delete;
+ /** Assignment operator (disabled) */
+ CbmRichDigitizer& operator=(const CbmRichDigitizer&) = delete;
- /**
- * \brief Inherited from FairTask.
+ /**
+ * \brief Detector system ID
+ * \return kRich
*/
+ ECbmModuleId GetSystemId() const { return ECbmModuleId::kRich; }
+
+ /** Inherited from FairTask */
virtual InitStatus Init();
- /**
- * \brief Inherited from FairTask.
- */
+ /** Inherited from FairTask */
virtual void Exec(Option_t* option);
- /**
- * \brief Inherited from FairTask.
- */
+ /** Inherited from FairTask */
virtual void Finish();
- /**
- * \brief Set crosstalk probability.
- */
- void SetCrossTalkProbability(Double_t crosstalk) { fCrossTalkProbability = crosstalk; }
-
- /**
- * \brief Set detector type
- */
+ /** Set detector type */
void SetDetectorType(CbmRichPmtTypeEnum detType) { fDetectorType = detType; }
- /**
- * \brief noise rate per McRichPoint / per pixel / per second :
- hofNoiseDigis = nofRichPoints * nofPixels * dT(50 ns) * (fNoiseDigiRate / 1.e9);
- */
- void SetNoiseDigiRate(Double_t noise) { fNoiseDigiRate = noise; }
+ /** Set time resolution for digis */
+ void SetTimeResolution(Double_t dt) { fTimeResolution = dt; }
/**
- * \brief Set collection efficiency for photoelectrons in PMT optics.
+ * \brief Set pixel dead time between signals (without time smearing).
+ * The resulting dead time at digi level is approximately:
+ * pixel dead time - 2*5*time resolution (5 sigma from gauss smearing for 2 digis each).
*/
- void SetCollectionEfficiency(Double_t collEff) { fPmt.SetCollectionEfficiency(collEff); }
+ void SetPixelDeadTime(Double_t dt) { fPixelDeadTime = dt; }
+
+ /** Set crosstalk probability */
+ void SetCrossTalkProbability(Double_t crosstalk) { fCrossTalkProbability = crosstalk; }
/**
- * \brief Set additional smearing of MC Points due to light scattering in mirror.
+ * \brief Set event noise rate per McRichPoint / per pixel / per second :
+ * nofNoiseSignals = nofRichPoints * nofPixels * event noise interval * (event noise rate / 1.e9);
*/
- //void SetSigmaMirror(Double_t sigMirror) {fSigmaMirror = sigMirror;}
+ void SetEventNoiseRate(Double_t noise) { fEventNoiseRate = noise; }
/**
- * \brief Set Time resolution.
+ * \brief Set event noise interval in ns.
+ * Add noise signals from: current event time to current event time + event noise interval.
*/
- void SetTimeResolution(Double_t dt) { fTimeResolution = dt; }
+ void SetEventNoiseInterval(Double_t dT) { fEventNoiseInterval = dT; }
+ /** Set dark rate per pixel in Hz */
+ void SetDarkRatePerPixel(Double_t darkRate) { fDarkRatePerPixel = darkRate; }
/**
- * \brief Set Pixel dead time.
+ * \brief Set charged particle cluster signal probability for direct neighbors.
+ * The probability for pixels further away decreases with 1/distance.
+ * Produced signals are limited to the same PMT as the source signal pixel and the defined cluster size.
*/
- void SetPixelDeadTime(Double_t dt) { fPixelDeadTime = dt; }
+ void SetClusterSignalProbability(Double_t intensity) { fClusterSignalProbability = intensity; }
/**
- * \brief Set Maximum nimber of hits per PMT cut.
+ * \brief Set cluster size for charged particle clusters.
+ * Resulting cluster size : (2*size + 1)*(2*size + 1) pixels.
*/
- void SetMaxNofHitsPerPmtCut(Double_t nofHits) { fMaxNofHitsPerPmtCut = nofHits; }
+ void SetClusterSize(UInt_t size) { fClusterSize = size; }
+
+ /** Set collection efficiency for photoelectrons in PMT optics */
+ void SetCollectionEfficiency(Double_t collEff) { fPmt.SetCollectionEfficiency(collEff); }
+
+ /** Set additional smearing of MC Points due to light scattering in mirror */
+ // void SetSigmaMirror(Double_t sigMirror) {fSigmaMirror = sigMirror;}
/**
* \brief Set if you want to shift z MC point value (workaround for GEANT4).
@@ -121,83 +121,82 @@ public:
private:
- Int_t fEventNum;
-
- TClonesArray* fRichPoints; // RICH MC points
- TClonesArray* fRichDigis; // RICH digis (output array)
- TClonesArray* fRichDigiMatches; // RICH digi matches (output array)
- TClonesArray* fMcTracks; // Monte-Carlo tracks
-
- Double_t fNofPoints; // total number of MCPoints processed
- Double_t fNofDigis; // total number of digis created
- Double_t fTimeTot; // sum of execution time
-
- CbmRichPmt fPmt;
- Double_t fCrossTalkProbability; // probability of the crosstalk for direct neighbor for one pixel
- Double_t fNoiseDigiRate; // noise rate per McRichPoint / per pixel / per second :
- // hofNoiseDigis = nofRichPoints * nofPixels * dT(50 ns) * (fNoiseDigiRate / 1.e9);
- CbmRichPmtTypeEnum fDetectorType;
- Int_t fMaxNofHitsPerPmtCut; // maximum number of hits which can be registered per PMT per event.
- // If more then the whole PMT is skipped
-
- map<Int_t, pair<CbmRichDigi*, CbmMatch*>> fDataMap;
-
- Double_t fTimeResolution; // in ns
- Double_t fDarkRatePerPixel; // dark rate per pixel in Hz
- Double_t fPixelDeadTime; // in ns, during this time pixel can not be fired
- map<Int_t, Double_t> fFiredPixelsMap; // first: pixel address, second: last fired time.
- Bool_t
- fDoZShift; // Set if you want to shift z MC point value (workaround for GEANT4). Must be set to true if one runs full RICH geoemtry with GEANT4, fot mCBM set to false
-
- /*
- * \brief Add crasstalk digis to the output array for the digi assuming fCrossTalkProbability
- */
- void AddCrossTalkDigis(CbmRichDigi* digi);
+ Int_t fEventNum {};
- /*
- * \brief Add fNofNoiseDigis number of digis.
- */
- void AddNoiseDigis(Int_t eventNum, Int_t inputNum);
+ TClonesArray* fRichPoints {nullptr}; // RICH MC points
+ TClonesArray* fRichDigis {nullptr}; // RICH digis (output array)
+ TClonesArray* fRichDigiMatches {nullptr}; // RICH digi matches (output array)
+ TClonesArray* fMcTracks {nullptr}; // Monte-Carlo tracks
- /*
- * \brief Process CbmRichPoint. Main method which is calle dfor all CbmRichPoints.
- */
- void ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t eventNum, Int_t inputNum);
+ Int_t fNofPoints {}; // total number of MCPoints processed
+ Int_t fNofDigis {}; // total number of digis created
+ Double_t fTimeTot {}; // sum of execution time
- /*
- * \brief Add all the fired digis to the output array
- * \@value Number of digis written
- */
- Int_t AddDigisToOutputArray();
+ CbmRichPmt fPmt {};
+ CbmRichPmtTypeEnum fDetectorType {CbmRichPmtTypeCosy17NoWls};
+
+ Double_t fTimeResolution {1.}; // in ns, time resolution for digis
+ Double_t fPixelDeadTime {50.}; // in ns, deadtime for one pixel
+ Double_t fCrossTalkProbability {0.02}; // probability of the crosstalk for direct neighbor for one pixel
+ Double_t fEventNoiseRate {5.}; // noise rate per McRichPoint / per pixel / per second
+ Double_t fEventNoiseInterval {50.}; // in ns, interval to generate event noise signals
+ Double_t fDarkRatePerPixel {1000.}; // dark rate per pixel in Hz
+ Double_t fClusterSignalProbability {0.37}; // cluster intensity
+ UInt_t fClusterSize {0}; // cluster size
+
+ std::map<Int_t, std::vector<std::pair<Double_t, CbmLink*>>>
+ fSignalBuffer {}; // first: pixel address, second: signal buffer
- /*
+ std::map<Int_t, Double_t> fPixelsLastFiredTime {}; // first: pixel address, second: last fired time
+ Bool_t fDoZShift {
+ true}; // Set if you want to shift z MC point value (workaround for GEANT4). Must be set to true if one runs full RICH geoemtry with GEANT4, fot mCBM set to false
+
+ /**
* \brief Process current MC event.
- * \value Number of processed RichPoints
+ * \return Number of processed RichPoints.
*/
Int_t ProcessMcEvent();
- /*
- * \brief Generate noise between events.
- */
- void GenerateNoiseBetweenEvents(Double_t oldEventTime, Double_t newEventTime);
+ /** Process CbmRichPoint. Main method which is called for all CbmRichPoints. */
+ void ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t eventNum, Int_t inputNum);
+
+ /** Add signal to signal buffer */
+ void AddSignalToBuffer(Int_t address, Double_t time, const CbmLink& link);
- /*
- *
+ /**
+ * \brief Add crosstalk assuming fCrossTalkProbability.
+ * Only add maximum one cross talk signal per MC point.
+ * Limited to the same PMT as the source MC point.
*/
- void AddDigi(Int_t address, Double_t time, const CbmLink& link);
+ void AddCrossTalk(Int_t address, Double_t time, const CbmLink& link);
+ /**
+ * \brief Add additional signals to nearby pixels if a charged particle
+ * directly passes through the PMT, given fClusterSignalProbability and fClusterSize.
+ * Limited to the same PMT. Currently independant of mass, momentum, etc.
+ */
+ void AddChargedParticleCluster(Int_t address, Double_t time, Int_t eventNum, Int_t inputNum);
/**
- * \brief Copy constructor.
+ * \brief Add additional noise for each event based on fEventNoiseRate
+ * and fEventNoiseInterval.
*/
- CbmRichDigitizer(const CbmRichDigitizer&);
+ void AddEventNoise(Int_t eventNum, Int_t inputNum);
+
+ /** Add noise between events based on fDarkRatePerPixel */
+ void AddDarkRateNoise(Double_t oldEventTime, Double_t newEventTime);
/**
- * \brief Assignment operator.
+ * \brief Process signals in all buffers until processTime.
+ * New Digis are only created until processTime - fPixelDeadTime,
+ * since potential links are not added to the buffer yet.
+ * Links are potentially added until processTime.
+ * \param processTime time in ns
+ * \return Number of Digis sent to DAQ
*/
- CbmRichDigitizer& operator=(const CbmRichDigitizer&);
+ Int_t ProcessBuffers(Double_t processTime);
- ClassDef(CbmRichDigitizer, 2)
+ ClassDef(CbmRichDigitizer, 3)
};
#endif
diff --git a/sim/response/base/CbmDigitization.cxx b/sim/response/base/CbmDigitization.cxx
index 5cbc82cacaa7ac2c8e787e3f1bafa77464256747..20590bab06fec3b2f4670569841a6d181f5567f1 100644
--- a/sim/response/base/CbmDigitization.cxx
+++ b/sim/response/base/CbmDigitization.cxx
@@ -129,7 +129,7 @@ Int_t CbmDigitization::CreateDefaultDigitizers()
{
std::cout << "Create default digitisers" << std::endl;
- stringstream ss;
+ std::stringstream ss;
ss << fName << ": Create default digitisers: ";
Int_t nDigis = 0;