diff --git a/sim/detectors/rich/CbmRich.cxx b/sim/detectors/rich/CbmRich.cxx
index 4777234decbb649a5b0a38796278b18f15b3988a..99d2d269a070ecd2cbfcd62e738b96b07dc301ac 100644
--- a/sim/detectors/rich/CbmRich.cxx
+++ b/sim/detectors/rich/CbmRich.cxx
@@ -1,44 +1,35 @@
-/* Copyright (C) 2006-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+/* Copyright (C) 2006-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: claudia Hoehne [committer], Andrey Lebedev, Florian Uhlig, Semen Lebedev */
+ Authors: claudia Hoehne [committer], Andrey Lebedev, Florian Uhlig, Semen Lebedev, Martin Beyer */
#include "CbmRich.h"
#include "CbmGeometryUtils.h"
+#include "CbmRichDigiMapManager.h"
#include "CbmRichPoint.h"
#include "CbmStack.h"
-#include "FairGeoInterface.h"
-#include "FairGeoLoader.h"
-#include "FairGeoNode.h"
-#include "FairGeoRootBuilder.h"
-#include "FairRootManager.h"
-#include "FairRun.h"
-#include "FairRuntimeDb.h"
-#include "FairVolume.h"
-
-//#include "Math/Interpolator.h"
-#include "FairGeoBuilder.h"
-#include "FairGeoMedia.h"
-
-#include "TClonesArray.h"
-#include "TGDMLParse.h"
-#include "TGeoMCGeometry.h"
-#include "TGeoManager.h"
-#include "TGeoMatrix.h"
-#include "TGeoMedium.h"
-#include "TGeoNode.h"
-#include "TLorentzVector.h"
-#include "TObjArray.h"
-#include "TParticle.h"
-#include "TParticlePDG.h"
-#include "TVirtualMC.h"
-#include <TFile.h>
+#include <FairGeoInterface.h>
+#include <FairGeoLoader.h>
+#include <FairRootManager.h>
+#include <FairRun.h>
+#include <FairVolume.h>
-#include <iostream>
+// #include <FairGeoBuilder.h>
+#include <FairGeoMedia.h>
-using std::cout;
-using std::endl;
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <TGDMLParse.h>
+#include <TGeoManager.h>
+#include <TGeoMatrix.h>
+#include <TGeoMedium.h>
+#include <TGeoNode.h>
+#include <TLorentzVector.h>
+#include <TObjArray.h>
+#include <TParticle.h>
+#include <TParticlePDG.h>
+#include <TVirtualMC.h>
std::map<TString, TGeoMedium*> CbmRich::fFixedMedia;
@@ -49,8 +40,8 @@ CbmRich::CbmRich()
, fRichPoints(new TClonesArray("CbmRichPoint"))
, fRichRefPlanePoints(new TClonesArray("CbmRichPoint"))
, fRichMirrorPoints(new TClonesArray("CbmRichPoint"))
- , fRotation(NULL)
- , fPositionRotation(NULL)
+ , fRotation(nullptr)
+ , fPositionRotation(nullptr)
{
fVerboseLevel = 1;
}
@@ -71,22 +62,26 @@ CbmRich::CbmRich(const char* name, Bool_t active, Double_t px, Double_t py, Doub
CbmRich::~CbmRich()
{
- if (NULL != fRichPoints) {
+ if (fRichPoints) {
fRichPoints->Delete();
delete fRichPoints;
}
- if (NULL != fRichRefPlanePoints) {
+ if (fRichRefPlanePoints) {
fRichRefPlanePoints->Delete();
delete fRichRefPlanePoints;
}
- if (NULL != fRichMirrorPoints) {
+ if (fRichMirrorPoints) {
fRichMirrorPoints->Delete();
delete fRichMirrorPoints;
}
}
-void CbmRich::Initialize() { FairDetector::Initialize(); }
+void CbmRich::Initialize()
+{
+ FairDetector::Initialize();
+ CbmRichDigiMapManager::GetInstance();
+}
Bool_t CbmRich::IsSensitive(const std::string& name)
{
@@ -107,9 +102,9 @@ Bool_t CbmRich::ProcessHits(FairVolume* vol)
Int_t iVol = vol->getMCid();
TString volName = TString(vol->GetName());
+ // Treat MAPMT pixels
if (volName.Contains("pmt_pixel")) {
if (gMC->IsTrackEntering()) {
- // cout << gGeoManager->GetPath() << endl;
TParticle* part = gMC->GetStack()->GetCurrentTrack();
Double_t charge = part->GetPDG()->Charge() / 3.;
Int_t trackID = gMC->GetStack()->GetCurrentTrackNumber();
@@ -117,30 +112,37 @@ Bool_t CbmRich::ProcessHits(FairVolume* vol)
Double_t length = gMC->TrackLength();
Double_t eLoss = gMC->Edep();
+ // For pmt pixels use unique pixel address as detId instead of volume id
+ const std::string path = std::string(gMC->CurrentVolPath());
+ Int_t pixelAddress = CbmRichDigiMapManager::GetInstance().GetPixelAddressByPath(path);
+ if (pixelAddress == -1) {
+ LOG(error) << "CbmRich::ProcessHits() Pixel address is not found for path " << path;
+ return kFALSE;
+ }
+
TLorentzVector tPos, tMom;
gMC->TrackPosition(tPos);
gMC->TrackMomentum(tMom);
if (pdgCode == 50000050) { // Cherenkovs only
-
- AddHit(trackID, iVol, TVector3(tPos.X(), tPos.Y(), tPos.Z()), TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time,
- length, eLoss);
+ AddHit(trackID, pixelAddress, TVector3(tPos.X(), tPos.Y(), tPos.Z()), TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),
+ time, length, eLoss);
// Increment number of RichPoints for this track
- CbmStack* stack = (CbmStack*) gMC->GetStack();
+ CbmStack* stack = static_cast<CbmStack*>(gMC->GetStack());
stack->AddPoint(ECbmModuleId::kRich);
return kTRUE;
}
else {
if (charge == 0.) {
- return kFALSE; // no neutrals
+ return kFALSE; // No neutrals
}
- else { // charged particles
- AddHit(trackID, iVol, TVector3(tPos.X(), tPos.Y(), tPos.Z()), TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time,
- length, eLoss);
+ else { // Charged particles
+ AddHit(trackID, pixelAddress, TVector3(tPos.X(), tPos.Y(), tPos.Z()),
+ TVector3(tMom.Px(), tMom.Py(), tMom.Pz()), time, length, eLoss);
// Increment number of RichPoints for this track
- CbmStack* stack = (CbmStack*) gMC->GetStack();
+ CbmStack* stack = static_cast<CbmStack*>(gMC->GetStack());
stack->AddPoint(ECbmModuleId::kRich);
return kTRUE;
}
@@ -148,10 +150,13 @@ Bool_t CbmRich::ProcessHits(FairVolume* vol)
}
}
- // Treat imaginary plane in front of the mirrors: Only charged particles at entrance
+ // Treat imaginary plane in front of the mirrors
+ // Two modes based on fRegisterPhotonsOnSensitivePlane:
+ // 1. [true] register only cherenkov photons
+ // Use case: check ingoing cherenkov photons into the mirror and outgoing cherenkov photons from the mirror
+ // 2. [false] register all charged particles
+ // Use case: compare track extrapolation with MC track
if (volName.Contains("sens_plane")) {
- // cout << volName << endl;
- // Collecting points of tracks and imaginary plane intersection
if (gMC->IsTrackEntering()) {
TParticle* part = gMC->GetStack()->GetCurrentTrack();
Double_t charge = part->GetPDG()->Charge() / 3.;
@@ -162,16 +167,16 @@ Bool_t CbmRich::ProcessHits(FairVolume* vol)
Double_t time = gMC->TrackTime() * 1.0e09;
Double_t length = gMC->TrackLength();
Double_t eLoss = gMC->Edep();
- TLorentzVector tPos, tMom;
+ TLorentzVector tPos, tMom;
gMC->TrackPosition(tPos);
gMC->TrackMomentum(tMom);
AddRefPlaneHit(trackID, iVol, TVector3(tPos.X(), tPos.Y(), tPos.Z()), TVector3(tMom.Px(), tMom.Py(), tMom.Pz()),
time, length, eLoss);
- //Increment number of RefPlanePoints for this track
- CbmStack* stack = (CbmStack*) gMC->GetStack();
+ // Increment number of RefPlanePoints for this track
+ CbmStack* stack = static_cast<CbmStack*>(gMC->GetStack());
stack->AddPoint(ECbmModuleId::kRef);
return kTRUE;
}
@@ -181,23 +186,18 @@ Bool_t CbmRich::ProcessHits(FairVolume* vol)
}
}
- // Treat mirror points
- Bool_t isMirror = (volName.Contains("mirror_tile_type"));
- if (isMirror) {
-
- // Collecting points
+ // Treat mirror
+ if (volName.Contains("mirror_tile_type")) {
if (gMC->IsTrackEntering()) {
TParticle* part = gMC->GetStack()->GetCurrentTrack();
Double_t charge = part->GetPDG()->Charge() / 3.;
if (charge != 0.) {
-
- Int_t trackID = gMC->GetStack()->GetCurrentTrackNumber();
-
+ Int_t trackID = gMC->GetStack()->GetCurrentTrackNumber();
Double_t time = gMC->TrackTime() * 1.0e09;
Double_t length = gMC->TrackLength();
Double_t eLoss = gMC->Edep();
- TLorentzVector tPos, tMom;
+ TLorentzVector tPos, tMom;
gMC->TrackPosition(tPos);
gMC->TrackMomentum(tMom);
@@ -229,7 +229,7 @@ TClonesArray* CbmRich::GetCollection(Int_t iColl) const
if (iColl == 0) return fRichPoints;
if (iColl == 1) return fRichRefPlanePoints;
if (iColl == 2) return fRichMirrorPoints;
- return NULL;
+ return nullptr;
}
void CbmRich::Print(Option_t*) const
@@ -254,10 +254,10 @@ void CbmRich::CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset)
{
Int_t nEntries = cl1->GetEntriesFast();
LOG(info) << "CbmRich: " << nEntries << " entries to add.";
- TClonesArray& clref = *cl2;
- CbmRichPoint* oldpoint = NULL;
+ TClonesArray& clref = *cl2;
+ CbmRichPoint* oldpoint {nullptr};
for (Int_t i = 0; i < nEntries; i++) {
- oldpoint = (CbmRichPoint*) cl1->At(i);
+ oldpoint = static_cast<CbmRichPoint*>(cl1->At(i));
Int_t index = oldpoint->GetTrackID() + offset;
oldpoint->SetTrackID(index);
new (clref[fPosIndex]) CbmRichPoint(*oldpoint);
@@ -319,9 +319,9 @@ void CbmRich::SetRichGlassPropertiesForGeant4()
2.6955, 2.6382, 2.5832, 2.5305, 2.4799, 2.4313, 2.3845, 2.3395, 2.2962, 2.2544, 2.2142, 2.1753, 2.1378,
2.1016, 2.0666, 2.0327, 1.9999, 1.9682, 1.9374, 1.9076, 1.8787, 1.8507, 1.8234, 1.797, 1.7713, 1.7464,
1.7221, 1.6985, 1.6756, 1.6533, 1.6315, 1.6103, 1.5897, 1.5695, 1.5499};
- // transform to GeV
- for (size_t i = 0; i < energyAr.size(); i++) {
- energyAr[i] = 1.e-9 * energyAr[i];
+ // Transform to GeV
+ for (Double_t& energy : energyAr) {
+ energy *= 1.e-9;
}
std::vector<Double_t> reflectivityAr = {
@@ -331,12 +331,13 @@ void CbmRich::SetRichGlassPropertiesForGeant4()
0.12021, 0.11909, 0.11783, 0.12257, 0.1215, 0.12199, 0.12494, 0.13101, 0.12089, 0.12284, 0.12569, 0.13136, 0.13307,
0.13705, 0.13844, 0.13753, 0.14416, 0.14761, 0.14953, 0.15218, 0.15315, 0.15719};
- for (size_t i = 0; i < reflectivityAr.size(); i++) {
- reflectivityAr[i] = 1. - reflectivityAr[i];
+ for (Double_t& reflectivity : reflectivityAr) {
+ reflectivity = 1. - reflectivity;
}
gMC->DefineOpSurface("RICHglassSurf", kGlisur, kDielectric_metal, kPolished, 0.0);
- gMC->SetMaterialProperty("RICHglassSurf", "REFLECTIVITY", energyAr.size(), energyAr.data(), reflectivityAr.data());
+ gMC->SetMaterialProperty("RICHglassSurf", "REFLECTIVITY", static_cast<Int_t>(energyAr.size()), energyAr.data(),
+ reflectivityAr.data());
for (int i = 0; i < 10; i++) {
if (gGeoManager->FindVolumeFast(("mirror_tile_type" + std::to_string(i)).c_str()) == nullptr) continue;
@@ -349,7 +350,7 @@ void CbmRich::ConstructGdmlGeometry(TGeoMatrix* geoMatrix)
{
TFile* old = gFile;
TGDMLParse parser;
- TGeoVolume* gdmlTop;
+ TGeoVolume* gdmlTop {nullptr};
// Before importing GDML
Int_t maxInd = gGeoManager->GetListOfMedia()->GetEntries() - 1;
@@ -361,10 +362,10 @@ void CbmRich::ConstructGdmlGeometry(TGeoMatrix* geoMatrix)
// gGeoManager->GetListOfMedia()->At(i)->Dump();
// After importing GDML
Int_t j = gGeoManager->GetListOfMedia()->GetEntries() - 1;
- Int_t curId;
- TGeoMedium* m;
+ Int_t curId {};
+ TGeoMedium* m {nullptr};
do {
- m = (TGeoMedium*) gGeoManager->GetListOfMedia()->At(j);
+ m = static_cast<TGeoMedium*>(gGeoManager->GetListOfMedia()->At(j));
curId = m->GetId();
m->SetId(curId + maxInd);
j--;
@@ -379,7 +380,7 @@ void CbmRich::ConstructGdmlGeometry(TGeoMatrix* geoMatrix)
ExpandNodeForGdml(gGeoManager->GetTopVolume()->GetNode(gGeoManager->GetTopVolume()->GetNdaughters() - 1));
for (Int_t k = maxInd + 1; k < newMaxInd + 1; k++) {
- TGeoMedium* medToDel = (TGeoMedium*) (gGeoManager->GetListOfMedia()->At(maxInd + 1));
+ TGeoMedium* medToDel = static_cast<TGeoMedium*>(gGeoManager->GetListOfMedia()->At(maxInd + 1));
LOG(debug) << " removing media " << medToDel->GetName() << " with id " << medToDel->GetId() << " (k=" << k
<< ")";
gGeoManager->GetListOfMedia()->Remove(medToDel);
@@ -429,7 +430,7 @@ void CbmRich::ExpandNodeForGdml(TGeoNode* node)
gGeoManager->RemoveMaterial(ind);
LOG(debug) << " removing material " << curMat->GetName() << " with index " << ind;
for (Int_t i = ind; i < gGeoManager->GetListOfMaterials()->GetEntries(); i++) {
- TGeoMaterial* m = (TGeoMaterial*) gGeoManager->GetListOfMaterials()->At(i);
+ TGeoMaterial* m = static_cast<TGeoMaterial*>(gGeoManager->GetListOfMaterials()->At(i));
m->SetIndex(m->GetIndex() - 1);
}
@@ -446,10 +447,10 @@ void CbmRich::ExpandNodeForGdml(TGeoNode* node)
FairGeoLoader* geoLoad = FairGeoLoader::Instance();
FairGeoInterface* geoFace = geoLoad->getGeoInterface();
FairGeoMedia* geoMediaBase = geoFace->getMedia();
- // FairGeoBuilder* geobuild = geoLoad->getGeoBuilder();
+ // FairGeoBuilder* geobuild = geoLoad->getGeoBuilder();
FairGeoMedium* curMedInGeo = geoMediaBase->getMedium(medName);
- if (curMedInGeo == 0) {
+ if (curMedInGeo == nullptr) {
LOG(fatal) << medName << " Media not found in Geo file.";
//! This should not happen.
//! This means that somebody uses material in GDML that is not in the media.geo file.
@@ -457,14 +458,14 @@ void CbmRich::ExpandNodeForGdml(TGeoNode* node)
}
else {
LOG(debug) << " Found media in Geo file" << medName;
- // Int_t nmed = geobuild->createMedium(curMedInGeo);
- fFixedMedia[medName] = (TGeoMedium*) gGeoManager->GetListOfMedia()->Last();
+ // Int_t nmed = geobuild->createMedium(curMedInGeo);
+ fFixedMedia[medName] = static_cast<TGeoMedium*>(gGeoManager->GetListOfMedia()->Last());
gGeoManager->RemoveMaterial(curMatOfMedInGeoManager->GetIndex());
LOG(debug) << " removing material " << curMatOfMedInGeoManager->GetName() << " with index "
<< curMatOfMedInGeoManager->GetIndex();
for (Int_t i = curMatOfMedInGeoManager->GetIndex(); i < gGeoManager->GetListOfMaterials()->GetEntries();
i++) {
- TGeoMaterial* m = (TGeoMaterial*) gGeoManager->GetListOfMaterials()->At(i);
+ TGeoMaterial* m = static_cast<TGeoMaterial*>(gGeoManager->GetListOfMaterials()->At(i));
m->SetIndex(m->GetIndex() - 1);
}
}
@@ -487,16 +488,16 @@ void CbmRich::ExpandNodeForGdml(TGeoNode* node)
curVol->SetMedium(fFixedMedia[medName]);
gGeoManager->SetAllIndex();
- // gGeoManager->GetListOfMaterials()->Print();
- // gGeoManager->GetListOfMedia()->Print();
+ // gGeoManager->GetListOfMaterials()->Print();
+ // gGeoManager->GetListOfMedia()->Print();
}
//! Recursevly go down the tree of nodes
if (curVol->GetNdaughters() != 0) {
TObjArray* NodeChildList = curVol->GetNodes();
- TGeoNode* curNodeChild;
+ TGeoNode* curNodeChild {nullptr};
for (Int_t j = 0; j < NodeChildList->GetEntriesFast(); j++) {
- curNodeChild = (TGeoNode*) NodeChildList->At(j);
+ curNodeChild = static_cast<TGeoNode*>(NodeChildList->At(j));
ExpandNodeForGdml(curNodeChild);
}
}
diff --git a/sim/detectors/rich/CbmRich.h b/sim/detectors/rich/CbmRich.h
index 773523db8f1361cff21072594f47b6f2d63ea2b2..98eb16afd71f9573eedf3ba757863fea178e8845 100644
--- a/sim/detectors/rich/CbmRich.h
+++ b/sim/detectors/rich/CbmRich.h
@@ -14,24 +14,21 @@
#ifndef CBM_RICH
#define CBM_RICH
+#include <FairDetector.h>
-#include "FairDetector.h"
-
-#include "Rtypes.h"
-#include "TGeoMatrix.h"
-#include "TString.h"
-#include "TVector3.h"
-#include <TClonesArray.h> // for ROOTCLING
+#include <Rtypes.h>
+#include <TClonesArray.h>
+#include <TString.h>
+#include <TVector3.h>
#include <map>
-class CbmRichRefPlanePoint;
class CbmRichPoint;
-class CbmRichMirrorPoint;
class FairVolume;
class TGeoMatrix;
class TGeoNode;
class TGeoMedium;
+
/**
* \class CbmRich
*
@@ -41,7 +38,6 @@ class TGeoMedium;
* \date 2004
**/
class CbmRich : public FairDetector {
-
public:
/**
* \brief Default constructor.
@@ -63,26 +59,23 @@ public:
Double_t pz = 258.75, // Z coordinate for v16a = 270, for v17a = 258.75, for v18a = 0
Double_t rx = 0., Double_t ry = 0., Double_t rz = 0.);
-
/**
* \brief Destructor.
*/
virtual ~CbmRich();
-
/**
* \brief Initialize detector. Stores volume IDs for RICH detector and mirror.
*/
virtual void Initialize();
-
/**
- * \brief Defines the action to be taken when a step is inside the
- * active volume. Creates CbmRichPoints and CbmRichMirrorPoints and adds
- * them to the collections.
+ * \brief Defines the action to be taken when a step is inside the active volume.
+ * Creates CbmRichPoints for MAPMT pixels, sensitive reference plane, mirror
+ * and adds them to seperate hit collections.
* \param[in] vol Pointer to the active volume.
*/
- virtual Bool_t ProcessHits(FairVolume* vol = 0);
+ virtual Bool_t ProcessHits(FairVolume* vol = nullptr);
/**
* \brief If verbosity level is set, print hit collection at the
@@ -90,31 +83,26 @@ public:
*/
virtual void EndOfEvent();
-
/**
* \brief Registers the hit collection in the ROOT manager.
*/
virtual void Register();
-
/**
* \brief Return hit collection.
*/
virtual TClonesArray* GetCollection(Int_t iColl) const;
-
/**
* \brief Screen output of hit collection.
*/
virtual void Print(Option_t*) const;
-
/**
* \brief Clears the hit collection.
*/
virtual void Reset();
-
/**
* \brief Copies the hit collection with a given track index offset.
* \param[in] cl1 Origin array.
@@ -123,14 +111,12 @@ public:
*/
virtual void CopyClones(TClonesArray* cl1, TClonesArray* cl2, Int_t offset);
-
/**
* \brief Construct geometry. Currently ROOT and ASCII formats are supported.
* The concrete method for geometry construction is called according to geometry file.
*/
virtual void ConstructGeometry();
-
/**
* \brief Construct geometry from GDML file.
* \param[in] geoMatrix Position and rotation of the RICH detector.
@@ -143,43 +129,42 @@ public:
*/
void ExpandNodeForGdml(TGeoNode* node);
-
/**
* \brief Put some optical properties.
*/
void ConstructOpGeometry();
/**
- * \brief Set Cherenkov propeties for RICH mirror.
- */
+ * \brief Set Cherenkov propeties for RICH mirror.
+ */
void SetRichGlassPropertiesForGeant4();
-
- /** Check whether a volume is sensitive.
- ** The decision is based on the volume name. Only used in case
- ** of ROOT geometry.
- ** @since 11.06.2012
- ** @param(name) Volume name
- ** @value kTRUE if volume is sensitive, else kFALSE
- **/
+ /**
+ * \brief Check whether a volume is sensitive.
+ * The decision is based on the volume name. Only used in case
+ * of ROOT geometry.
+ * \since 11.06.2012
+ * \param name Volume name
+ * \return kTRUE if volume is sensitive, else kFALSE
+ */
virtual Bool_t CheckIfSensitive(std::string name);
virtual Bool_t IsSensitive(const std::string& name);
- /*
+ /**
* \brief set fRegisterPhotonsOnSensitivePlane parameter
*/
void SetRegisterPhotonsOnSensitivePlane(Bool_t b) { fRegisterPhotonsOnSensitivePlane = b; }
private:
Int_t fPosIndex;
- // set to true if you want to register photons onto the sensitive gas plane,
+
+ // Set to true if you want to register photons onto the sensitive gas plane,
// if false then only charged particles are registered
Bool_t fRegisterPhotonsOnSensitivePlane;
TClonesArray* fRichPoints; // MC points onto the photodetector plane
- TClonesArray* fRichRefPlanePoints; // points on the reference plane
- TClonesArray* fRichMirrorPoints; // mirror points
-
+ TClonesArray* fRichRefPlanePoints; // Points on the reference plane
+ TClonesArray* fRichMirrorPoints; // Mirror points
// GDML geometry
static std::map<TString, TGeoMedium*> fFixedMedia; // List of media "repaired" after importing GMDL
diff --git a/sim/detectors/rich/CbmRichDigitizer.cxx b/sim/detectors/rich/CbmRichDigitizer.cxx
index dedeebf2c3e0175cf588d975c4994417ff097ecb..db6fe7fcc0a1ba8cbaa7258e70ae7ad2616f5dfb 100644
--- a/sim/detectors/rich/CbmRichDigitizer.cxx
+++ b/sim/detectors/rich/CbmRichDigitizer.cxx
@@ -20,18 +20,15 @@
#include "CbmRichGeoManager.h"
#include "CbmRichPoint.h"
-#include <FairEventHeader.h>
-#include <FairRunAna.h>
-#include <FairRunSim.h>
+#include <FairRootManager.h>
#include <Logger.h>
#include <TClonesArray.h>
-#include <TGeoManager.h>
#include <TMath.h>
#include <TRandom.h>
#include <TStopwatch.h>
-#include <cassert>
+#include <cstdlib>
#include <iomanip>
#include <iostream>
@@ -135,7 +132,7 @@ void CbmRichDigitizer::Finish()
LOG(info) << "=====================================";
LOG(info) << GetName() << ": Run summary";
LOG(info) << "Events processed : " << fEventNum;
- LOG(info) << "RichPoint / event : " << setprecision(1) << fNofPoints / Double_t(fEventNum);
+ LOG(info) << "RichPoint / event : " << setprecision(6) << 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";
@@ -150,7 +147,7 @@ Int_t CbmRichDigitizer::ProcessMcEvent()
<< " EventTime:" << fCurrentEventTime << " nofRichPoints:" << nofRichPoints;
for (Int_t j = 0; j < nofRichPoints; j++) {
- CbmRichPoint* point = (CbmRichPoint*) fRichPoints->At(j);
+ CbmRichPoint* point = static_cast<CbmRichPoint*>(fRichPoints->At(j));
ProcessPoint(point, j, fCurrentMCEntry, fCurrentInput);
}
@@ -161,44 +158,37 @@ Int_t CbmRichDigitizer::ProcessMcEvent()
void CbmRichDigitizer::ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t eventNum, Int_t inputNum)
{
- // LOG(info) << "ProcessPoint " << pointId;
- // TGeoNode* node = gGeoManager->FindNode(point->GetX(), point->GetY(), point->GetZ());
- // workaround for GEANT4, probably boundary problem
- Double_t zNew = point->GetZ();
- if (fDoZShift) zNew = zNew - 0.001;
- gGeoManager->FindNode(point->GetX(), point->GetY(), zNew);
- std::string path(gGeoManager->GetPath());
- Int_t address = CbmRichDigiMapManager::GetInstance().GetPixelAddressByPath(path);
- if (address < 0) return;
+ Int_t address = point->GetDetectorID();
+ if (address == -1) {
+ LOG(error) << "CbmRichDigitizer::ProcessPoint: address == -1";
+ 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);
+ CbmMCTrack* mcTrack = static_cast<CbmMCTrack*>(fMcTracks->At(trackId));
+ if (!mcTrack) return;
+ Int_t gcode = TMath::Abs(mcTrack->GetPdgCode());
Bool_t isDetected = false;
- // for photons weight with efficiency of PMT
+ // For photons weight with quantum efficiency of PMT
if (gcode == 50000050) {
TVector3 mom;
point->Momentum(mom);
- Double_t momTotal = sqrt(mom.Px() * mom.Px() + mom.Py() * mom.Py() + mom.Pz() * mom.Pz());
+ Double_t momTotal = TMath::Sqrt(mom.Px() * mom.Px() + mom.Py() * mom.Py() + mom.Pz() * mom.Pz());
isDetected = fPmt.isPhotonDetected(fDetectorType, momTotal);
}
- else { // if not photon
- // worst case: assume that all charged particles crossing
+ else { // If not photon
+ // Worst case: assume that all charged particles crossing
// the PMTplane leave Cherenkov light in the PMTglass which will be detected
isDetected = true;
}
- Double_t time = fCurrentEventTime + point->GetTime();
if (isDetected) {
+ Double_t time = fCurrentEventTime + point->GetTime();
+ CbmLink link(1., pointId, eventNum, inputNum);
AddSignalToBuffer(address, time, link);
if (gcode == 50000050 && address > 0) AddCrossTalk(address, time, link);
- if (charge != 0. && address > 0) AddChargedParticleCluster(address, time, eventNum, inputNum);
+ if (mcTrack->GetCharge() != 0. && address > 0) AddChargedParticleCluster(address, time, eventNum, inputNum);
}
}
@@ -214,10 +204,11 @@ void CbmRichDigitizer::AddCrossTalk(Int_t address, Double_t time, const CbmLink&
/* 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.);
+ Double_t thresholdDirectDiagonal =
+ static_cast<Double_t>(directPixels.size())
+ / (static_cast<Double_t>(directPixels.size()) + static_cast<Double_t>(diagonalPixels.size()) / 4.);
if (gRandom->Uniform() < thresholdDirectDiagonal) {
crosstalkAddress = directPixels[gRandom->Integer(directPixels.size())];
} else {
@@ -230,23 +221,24 @@ void CbmRichDigitizer::AddCrossTalk(Int_t address, Double_t time, const CbmLink&
void CbmRichDigitizer::AddChargedParticleCluster(Int_t address, Double_t time, Int_t eventNum, Int_t inputNum)
{
- Int_t sourcePixelId = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(address)->fPixelId;
+ // pixelId % 8 is x index, pixelId / 8 is y index
+ auto sourcePixelIdx = std::div(CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(address)->fPixelId, 8);
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));
+ for (const Int_t& addr : neighbourAddresses) {
+ auto iPixelIdx = std::div(CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(addr)->fPixelId, 8);
+ Double_t d = TMath::Sqrt((sourcePixelIdx.rem - iPixelIdx.rem) * (sourcePixelIdx.rem - iPixelIdx.rem)
+ + (sourcePixelIdx.quot - iPixelIdx.quot) * (sourcePixelIdx.quot - iPixelIdx.quot));
if (gRandom->Uniform(1.) < fClusterSignalProbability / d) {
CbmLink link(1., -1, eventNum, inputNum);
- AddSignalToBuffer(neighbourAddresses[i], time, link);
+ AddSignalToBuffer(addr, time, link);
}
}
}
void CbmRichDigitizer::AddEventNoise(Int_t eventNum, Int_t inputNum)
{
- Int_t nofPixels = CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size();
+ Int_t nofPixels = static_cast<Int_t>(CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size());
Double_t nofRichPoints = fRichPoints->GetEntriesFast();
Double_t nofNoiseDigis = nofRichPoints * nofPixels * fEventNoiseInterval * (fEventNoiseRate / 1.e9);
@@ -262,7 +254,7 @@ void CbmRichDigitizer::AddEventNoise(Int_t eventNum, Int_t inputNum)
void CbmRichDigitizer::AddDarkRateNoise(Double_t oldEventTime, Double_t newEventTime)
{
- Int_t nofPixels = CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size();
+ Int_t nofPixels = static_cast<Int_t>(CbmRichDigiMapManager::GetInstance().GetPixelAddresses().size());
Double_t dT = newEventTime - oldEventTime;
Double_t nofNoisePixels = nofPixels * dT * (fDarkRatePerPixel / 1.e9);
@@ -314,8 +306,10 @@ Int_t CbmRichDigitizer::ProcessBuffers(Double_t processTime)
digi = new CbmRichDigi();
digi->SetAddress(dm.first);
digi->SetTime(digiTime);
- digiMatch = new CbmMatch();
- digiMatch->AddLink(*link);
+ if (fCreateMatches) { // If fCreateMatches is false, SendData does not take over ownership via unique_ptr
+ digiMatch = new CbmMatch();
+ digiMatch->AddLink(*link);
+ }
fPixelsLastFiredTime[dm.first] = signalTime; // Deadtime is from signal to signal
}
else {
diff --git a/sim/detectors/rich/CbmRichDigitizer.h b/sim/detectors/rich/CbmRichDigitizer.h
index 6283766327bbe9157ecaa1ca91cf6d9534cad87b..9d4af0a302ed2c9e3788664b494fb713c66dad9e 100644
--- a/sim/detectors/rich/CbmRichDigitizer.h
+++ b/sim/detectors/rich/CbmRichDigitizer.h
@@ -26,7 +26,6 @@ class TClonesArray;
class CbmRichPoint;
class CbmLink;
-
/**
* \class CbmRichDigitizer
*
@@ -114,12 +113,6 @@ public:
/** 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).
- */
- void SetDoZShift(Bool_t doZShift) { fDoZShift = doZShift; }
-
-
private:
Int_t fEventNum {};
@@ -148,8 +141,6 @@ private:
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.