Skip to content
Snippets Groups Projects
Commit 4b459cc0 authored by Lukas Chlad's avatar Lukas Chlad Committed by Volker Friese
Browse files

FSD integration into SIM

parent 8243afa7
No related branches found
No related tags found
1 merge request!1262FSD integration into SIM
Pipeline #23733 passed
Showing
with 1004 additions and 217 deletions
......@@ -103,7 +103,8 @@
#pragma link C++ class CbmFsdDigiData + ;
#pragma link C++ class CbmFsdHit + ;
#pragma link C++ class CbmFsdPoint + ;
#pragma link C++ class CbmFsdAddress + ;
#pragma link C++ namespace CbmFsdAddress;
#pragma link C++ enum CbmFsdAddress::Level;
// --- data/global
#pragma link C++ class CbmGlobalTrack + ;
......
......
......@@ -14,12 +14,16 @@
#include <Logger.h> // for Logger, LOG
// ----- Default constructor -------------------------------------------
CbmFsdHit::CbmFsdHit() : CbmPixelHit(), fModuleID(-1), fEdep(-1)
CbmFsdHit::CbmFsdHit() : CbmPixelHit(), fUnitId(-1), fModuleId(-1), fEdep(-1)
{
SetType(kFSDHIT);
SetTime(0.);
}
CbmFsdHit::CbmFsdHit(int32_t module, double edep) : CbmPixelHit(), fModuleID(module), fEdep(edep)
CbmFsdHit::CbmFsdHit(int32_t unit, int32_t module, double edep)
: CbmPixelHit()
, fUnitId(unit)
, fModuleId(module)
, fEdep(edep)
{
SetType(kFSDHIT);
SetTime(0.);
......@@ -35,7 +39,8 @@ void CbmFsdHit::Print(Option_t*) const { LOG(info) << ToString(); }
std::string CbmFsdHit::ToString() const
{
std::stringstream ss;
ss << "module : " << fModuleID << "position: [" << GetX() << "," << GetY() << "," << GetZ() << "] "
ss << "unit : " << fUnitId << "module : " << fModuleId << "position: [" << GetX() << "," << GetY() << "," << GetZ()
<< "] "
<< " ELoss " << fEdep;
return ss.str();
}
......
......
......@@ -28,7 +28,7 @@ public:
/** Default constructor **/
CbmFsdHit();
CbmFsdHit(int32_t module, double edep);
CbmFsdHit(int32_t unit, int32_t module, double edep);
/** Destructor **/
......@@ -40,8 +40,11 @@ public:
double GetEdep() const { return fEdep; }
void SetEdep(double edep) { fEdep = edep; }
int32_t GetModuleID() const { return fModuleID; }
void SetModuleID(int32_t mod) { fModuleID = mod; }
int32_t GetModuleId() const { return fModuleId; }
void SetModuleId(int32_t mod) { fModuleId = mod; }
int32_t GetUnitId() const { return fUnitId; }
void SetUnitId(int32_t unit) { fUnitId = unit; }
void Print(Option_t* = "") const;
......@@ -50,7 +53,8 @@ public:
private:
/** Data members **/
int32_t fModuleID;
int32_t fModuleId;
int32_t fUnitId;
double fEdep;
......
......
......@@ -10,7 +10,7 @@
#include <sstream> // for stringstream
// ----- Default constructor -------------------------------------------
CbmFsdPoint::CbmFsdPoint() : FairMCPoint(), fModuleID(0) {}
CbmFsdPoint::CbmFsdPoint() : FairMCPoint() {}
// -------------------------------------------------------------------------
......@@ -18,7 +18,6 @@ CbmFsdPoint::CbmFsdPoint() : FairMCPoint(), fModuleID(0) {}
CbmFsdPoint::CbmFsdPoint(int32_t trackID, int32_t detID, TVector3 pos, TVector3 mom, double tof, double length,
double eLoss)
: FairMCPoint(trackID, detID, pos, mom, tof, length, eLoss)
, fModuleID(0)
{
}
// -------------------------------------------------------------------------
......
......
......@@ -30,7 +30,7 @@ public:
/** Constructor with arguments
*@param trackID Index of MCTrack
*@param detID Detector ID
*@param pos Ccoordinates at entrance to active volume [cm]
*@param pos Coordinates at entrance to active volume [cm]
*@param mom Momentum of track at entrance [GeV]
*@param tof Time since event start [ns]
*@param length Track length since creation [cm]
......@@ -39,10 +39,6 @@ public:
CbmFsdPoint(int32_t trackID, int32_t detID, TVector3 pos, TVector3 mom, double tof, double length, double eLoss);
/** Copy constructor **/
// CbmFsdPoint(const CbmFsdPoint& point) { *this = point; };
/** Destructor **/
virtual ~CbmFsdPoint();
......@@ -50,17 +46,8 @@ public:
/** Output to screen **/
virtual void Print(const Option_t* opt) const;
/** Modifiers **/
void SetModuleID(int32_t mod) { fModuleID = mod; }
/** Accessors **/
int32_t GetModuleID() const { return fModuleID; }
std::string ToString() const;
private:
int32_t fModuleID; //number of module
ClassDef(CbmFsdPoint, 1)
};
......
......
......@@ -4,24 +4,27 @@ set(INCLUDE_DIRECTORIES
set(SRCS
CbmFsdContFact.cxx
CbmFsdDetectorMapManager.cxx
CbmFsdGeoHandler.cxx
CbmFsdDigiPar.cxx
)
set(LIBRARY_NAME CbmFsdBase)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
ROOT::Core
FairRoot::ParBase
)
set(PRIVATE_DEPENDENCIES
CbmData
FairLogger::FairLogger
ROOT::Core
ROOT::Geom
${VMCLIB}
)
generate_cbm_library()
# Install file which has no corresponding source file
install(FILES
CbmFsdDetectorData.h
DESTINATION include
set(HEADERS
CbmFsdDetectorSpecs.h
)
generate_cbm_library()
......@@ -11,6 +11,7 @@
#pragma link off all functions;
#pragma link C++ class CbmFsdContFact + ;
#pragma link C++ class CbmFsdDetectorMapManager + ;
#pragma link C++ class CbmFsdGeoHandler + ;
#pragma link C++ class CbmFsdDigiPar + ;
#endif
......@@ -13,6 +13,8 @@
*/
#include "CbmFsdContFact.h"
#include "CbmFsdDigiPar.h"
#include <FairContFact.h> // for FairContainer
#include <FairRuntimeDb.h> // for FairRuntimeDb
#include <Logger.h> // for LOG
......@@ -39,6 +41,13 @@ void CbmFsdContFact::setAllContainers()
{
/** Creates the Container objects with all accepted contexts and adds them to
* the list of containers for the FsdBase library.*/
FairContainer* p1 =
new FairContainer("CbmFsdDigiPar", "Digitization parameters for the FSD detector",
"Needed parameters to adjust FsdDigitizer according to the geometry and read-out propetries");
p1->addContext("Needed parameters to adjust FsdDigitizer according to the geometry and read-out propetries");
containers->Add(p1);
}
FairParSet* CbmFsdContFact::createContainer(FairContainer* c)
......@@ -48,7 +57,11 @@ FairParSet* CbmFsdContFact::createContainer(FairContainer* c)
* of this container, the name is concatinated with the context. */
const char* name = c->GetName();
LOG(info) << " -I container name " << name;
FairParSet* p = 0;
FairParSet* p = nullptr;
if (strcmp(name, "CbmFsdDigiPar") == 0) {
p = new CbmFsdDigiPar(c->getConcatName().Data(), c->GetTitle(), c->getContext());
}
return p;
}
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Lukas Chlad [committer] */
#include "CbmFsdDetectorMapManager.h"
#include "CbmFsdDetectorData.h" // for CbmFsdModuleData
#include <Logger.h> // for LOG, Logger
#include <TGeoBBox.h> // for TGeoBBox
#include <TGeoManager.h> // for TGeoManager, gGeoManager
#include <TGeoMatrix.h> // for TGeoMatrix
#include <TGeoNode.h> // for TGeoIterator, TGeoNode
#include <TGeoVolume.h> // for TGeoVolume
#include <TRandom.h> // for TRandom, gRandom
#include <string> // for operator<, stoul
#include <utility> // for pair
#include <stddef.h> // for size_t
using namespace std;
CbmFsdDetectorMapManager::CbmFsdDetectorMapManager() : fModulePathToIdMap(), fModuleIdToDataMap(), fModuleIds()
{
Init();
}
CbmFsdDetectorMapManager::~CbmFsdDetectorMapManager() {}
void CbmFsdDetectorMapManager::Init()
{
fModulePathToIdMap.clear();
fModuleIdToDataMap.clear();
fModuleIds.clear();
Int_t currentModuleId = 0;
TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
TGeoNode* curNode;
TString branchStr("fsd_"); // stored in path of nodes, fsd_geoTag_0
TString moduleStr("fsdmodule");
TString scintStr("scintillator");
geoIterator.Reset(); // safety to reset to "cave" befor the loop starts
while ((curNode = geoIterator())) {
TString nodePath;
geoIterator.GetPath(nodePath);
if (!nodePath.Contains(branchStr)) {
geoIterator.Skip(); // skip current branch when it is not FSD => should speed up
continue; // don't do anything for this branch
}
TString nodeName(curNode->GetName());
// FSD detector -> all modules names contain "fsdmodule"
if (TString(curNode->GetVolume()->GetName()).Contains(moduleStr)) {
const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
const Double_t* curNodeTr = curMatrix->GetTranslation();
string path = string(nodePath.Data());
TGeoVolume* curScintillator = nullptr;
for (int idn = 0; idn < curNode->GetNdaughters(); idn++) {
TGeoNode* daughterNode = curNode->GetDaughter(idn);
if (TString(daughterNode->GetVolume()->GetName()).Contains(scintStr)) {
curScintillator = daughterNode->GetVolume();
break;
}
}
const TGeoBBox* shape = (const TGeoBBox*) (curScintillator->GetShape());
fModulePathToIdMap.insert(pair<string, Int_t>(path, currentModuleId));
CbmFsdModuleData* moduleData = new CbmFsdModuleData();
moduleData->fX = curNodeTr[0];
moduleData->fY = curNodeTr[1];
moduleData->fZ = curNodeTr[2];
moduleData->dX = shape->GetDX();
moduleData->dY = shape->GetDY();
moduleData->dZ = shape->GetDZ();
moduleData->fId = currentModuleId;
fModuleIdToDataMap.insert(pair<Int_t, CbmFsdModuleData*>(moduleData->fId, moduleData));
fModuleIds.push_back(currentModuleId);
currentModuleId++;
}
}
LOG(info) << "CbmFsdDetectorMapManager is initialized";
LOG(info) << "fModulePathToIdMap.size() = " << fModulePathToIdMap.size();
LOG(info) << "fModuleIdToDataMap.size() = " << fModuleIdToDataMap.size();
}
Int_t CbmFsdDetectorMapManager::GetModuleIdByPath(const string& path)
{
std::map<string, Int_t>::iterator it;
it = fModulePathToIdMap.find(path);
if (it == fModulePathToIdMap.end()) return -1;
return it->second;
}
CbmFsdModuleData* CbmFsdDetectorMapManager::GetModuleDataById(Int_t modId)
{
std::map<Int_t, CbmFsdModuleData*>::iterator it;
it = fModuleIdToDataMap.find(modId);
if (it == fModuleIdToDataMap.end()) return nullptr;
return it->second;
}
vector<Int_t> CbmFsdDetectorMapManager::GetModuleIds() { return fModuleIds; }
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Lukas Chlad [committer] */
#ifndef CBMFSDDETECTORMAPMANAGER_H
#define CBMFSDDETECTORMAPMANAGER_H
#include <RtypesCore.h> // for Int_t, Double_t
#include <iostream> // for string
#include <map> // for map
#include <vector> // for vector
struct CbmFsdModuleData;
class CbmFsdDetectorMapManager {
private:
/** Default constructor **/
CbmFsdDetectorMapManager();
/** Default destructor **/
~CbmFsdDetectorMapManager();
public:
CbmFsdDetectorMapManager(const CbmFsdDetectorMapManager&) = delete;
CbmFsdDetectorMapManager& operator=(const CbmFsdDetectorMapManager&) = delete;
/**
* Return Instance of CbmFsdDetectorMapManager.
*/
static CbmFsdDetectorMapManager& GetInstance()
{
static CbmFsdDetectorMapManager fInstance;
return fInstance;
}
/*
* \brief Return digi moduleId by path to node.
*/
Int_t GetModuleIdByPath(const std::string& path);
/*
* \brief Return CbmFsdModuleData by digi moduleId.
*/
CbmFsdModuleData* GetModuleDataById(Int_t id);
/*
* \brief Return ids of all modules
*/
std::vector<Int_t> GetModuleIds();
private:
std::map<std::string, Int_t> fModulePathToIdMap;
std::map<Int_t, CbmFsdModuleData*> fModuleIdToDataMap;
std::vector<Int_t> fModuleIds; // vector of all module ids
/*
* \brief Initialize maps.
*/
void Init();
};
#endif /* CBMFSDDETECTORMAPMANAGER_H */
......@@ -2,14 +2,16 @@
SPDX-License-Identifier: GPL-3.0-only
Authors: Lukas Chlad [committer] */
#ifndef CBMFSDDETECTORDATA_H
#define CBMFSDDETECTORDATA_H
#ifndef CBMFSDDETECTORSPECS_H
#define CBMFSDDETECTORSPECS_H
#include <RtypesCore.h> // for Double_t, Int_t
#include <TString.h> // for TString
struct CbmFsdModuleData {
// module index
Int_t fId;
struct CbmFsdModuleSpecs {
// location indices
Int_t fUnitId;
Int_t fModuleId;
// module center position
Double_t fX;
Double_t fY;
......@@ -20,4 +22,13 @@ struct CbmFsdModuleData {
Double_t dZ;
};
#endif /* CBMFSDDETECTORDATA_H */
struct CbmFsdUnitSpecs {
// name of unit
TString fUnitName;
// unit id
Int_t fUnitId;
// number of modules in unit
Int_t fNumModules;
};
#endif /* CBMFSDDETECTORSPECS_H */
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Florian Uhlig, Lukas Chlad [committer] */
#include "CbmFsdDigiPar.h"
#include <FairParGenericSet.h> // for FairParGenericSet
#include <FairParamList.h> // for FairParamList
#include <Logger.h> // for LOG, Logger
ClassImp(CbmFsdDigiPar)
CbmFsdDigiPar::CbmFsdDigiPar(const char* name, const char* title, const char* context)
: FairParGenericSet(name, title, context)
, fNumPhotoDets(-1)
, fTimeResolution(-1.)
, fEnergyResolution(-1.)
{
detName = "Fsd";
}
CbmFsdDigiPar::~CbmFsdDigiPar(void) { clear(); }
void CbmFsdDigiPar::clear(void)
{
status = kFALSE;
resetInputVersions();
}
void CbmFsdDigiPar::putParams(FairParamList* l)
{
if (!l) { return; }
l->add("NumPhotoDets", fNumPhotoDets);
l->add("NumUnits", fNumUnits);
l->add("TimeResolution", fTimeResolution);
l->add("EnergyResolution", fEnergyResolution);
l->add("DeadTime", fDeadTime);
}
Bool_t CbmFsdDigiPar::getParams(FairParamList* l)
{
if (!l) { return kFALSE; }
LOG(debug2) << "Get the FSD digitization parameters.";
if (!l->fill("NumPhotoDets", &fNumPhotoDets)) return kFALSE;
if (!l->fill("NumUnits", &fNumUnits)) return kFALSE;
fTimeResolution.Set(fNumUnits);
fEnergyResolution.Set(fNumUnits);
fDeadTime.Set(fNumUnits);
if (!l->fill("TimeResolution", &fTimeResolution)) return kFALSE;
if (!l->fill("EnergyResolution", &fEnergyResolution)) return kFALSE;
if (!l->fill("DeadTime", &fDeadTime)) return kFALSE;
return kTRUE;
}
Double_t CbmFsdDigiPar::GetTimeResolution(Int_t iUnitId) const
{
if (iUnitId < fNumUnits) return fTimeResolution[iUnitId];
else
return -1.;
}
Double_t CbmFsdDigiPar::GetEnergyResolution(Int_t iUnitId) const
{
if (iUnitId < fNumUnits) return fEnergyResolution[iUnitId];
else
return -1.;
}
Double_t CbmFsdDigiPar::GetDeadTime(Int_t iUnitId) const
{
if (iUnitId < fNumUnits) return fDeadTime[iUnitId];
else
return -1.;
}
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Florian Uhlig, Lukas Chlad [committer] */
#ifndef CBMFSDDIGIPAR_H
#define CBMFSDDIGIPAR_H
#include <FairParGenericSet.h> // for FairParGenericSet
#include <Rtypes.h> // for THashConsistencyHolder, ClassDef
#include <RtypesCore.h> // for Int_t, Double_t
#include <TArrayD.h> // for TArrayD
class FairParamList;
class CbmFsdDigiPar : public FairParGenericSet {
public:
CbmFsdDigiPar(
const char* name = "CbmFsdDigiPar", const char* title = "Digitization parameters for the FSD detector",
const char* context = "Needed parameters to adjust FsdDigitizer according to the geometry and read-out propetries");
CbmFsdDigiPar(const CbmFsdDigiPar&) = delete;
CbmFsdDigiPar& operator=(const CbmFsdDigiPar&) = delete;
~CbmFsdDigiPar(void);
void clear(void);
void putParams(FairParamList*);
Bool_t getParams(FairParamList*);
Int_t GetNumPhotoDets() const { return fNumPhotoDets; }
Int_t GetNumUnits() const { return fNumUnits; }
Double_t GetTimeResolution(Int_t iUnitId) const;
Double_t GetEnergyResolution(Int_t iUnitId) const;
Double_t GetDeadTime(Int_t iUnitId) const;
private:
Int_t fNumPhotoDets; // number of photo detectors per module
Int_t fNumUnits; // number of units within given FSD geo version
TArrayD fTimeResolution; // value to smear the timing via gaussian
TArrayD fEnergyResolution; // value to smear the energy measurement via gaussian
TArrayD fDeadTime; // value to separate digis in time-based
ClassDef(CbmFsdDigiPar, 1)
};
#endif
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Lukas Chlad [committer] */
#include "CbmFsdGeoHandler.h"
#include "CbmFsdAddress.h" // for CbmFsdAddress
#include "CbmFsdDetectorSpecs.h" // for CbmFsdModuleSpecs
#include <Logger.h> // for LOG, Logger
#include <TGeoBBox.h> // for TGeoBBox
#include <TGeoManager.h> // for TGeoManager, gGeoManager
#include <TGeoMatrix.h> // for TGeoMatrix
#include <TGeoNode.h> // for TGeoIterator, TGeoNode
#include <TGeoVolume.h> // for TGeoVolume
#include <TVirtualMC.h> // for TVirtualMC, gMC
#include <string> // for operator<, stoul
#include <utility> // for pair
#include <stddef.h> // for size_t
using namespace std;
// constructor with initialization of maps
CbmFsdGeoHandler::CbmFsdGeoHandler() { InitMaps(); }
void CbmFsdGeoHandler::InitMaps()
{
fUnitIdToSpecsMap.clear();
fModuleIdToSpecsMap.clear();
Int_t currentModuleId = -1;
Int_t currentUnitId = -1;
TGeoIterator geoIterator(gGeoManager->GetTopNode()->GetVolume());
TGeoNode* curNode;
geoIterator.Reset(); // safety to reset to "cave" befor the loop starts
while ((curNode = geoIterator())) {
TString nodePath;
geoIterator.GetPath(nodePath);
if (!nodePath.Contains(fBranchStr)) {
geoIterator.Skip(); // skip current branch when it is not FSD => should speed up
continue; // don't do anything for this branch
}
TString nodeName(curNode->GetName());
if (nodeName.Contains(fUnitStr)) {
currentUnitId = curNode->GetNumber();
CbmFsdUnitSpecs* unitSpecs = new CbmFsdUnitSpecs();
unitSpecs->fUnitId = currentUnitId;
unitSpecs->fUnitName = static_cast<TString>(curNode->GetVolume()->GetName());
unitSpecs->fNumModules = curNode->GetNdaughters();
fUnitIdToSpecsMap.insert(pair<Int_t, CbmFsdUnitSpecs*>(currentUnitId, unitSpecs));
}
if (nodeName.Contains(fModuleStr)) {
currentModuleId = curNode->GetNumber();
const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
const Double_t* curNodeTr = curMatrix->GetTranslation();
TGeoVolume* curScintillator = nullptr;
for (int idn = 0; idn < curNode->GetNdaughters(); idn++) {
TGeoNode* daughterNode = curNode->GetDaughter(idn);
if (TString(daughterNode->GetName()).Contains(fActiveMatStr)) {
curScintillator = daughterNode->GetVolume();
break;
}
}
const TGeoBBox* shape = (const TGeoBBox*) (curScintillator->GetShape());
CbmFsdModuleSpecs* moduleSpecs = new CbmFsdModuleSpecs();
moduleSpecs->fX = curNodeTr[0];
moduleSpecs->fY = curNodeTr[1];
moduleSpecs->fZ = curNodeTr[2];
moduleSpecs->dX = shape->GetDX();
moduleSpecs->dY = shape->GetDY();
moduleSpecs->dZ = shape->GetDZ();
moduleSpecs->fModuleId = currentModuleId;
moduleSpecs->fUnitId = currentUnitId;
fModuleIdToSpecsMap.insert(pair<Int_t, CbmFsdModuleSpecs*>(currentModuleId, moduleSpecs));
}
}
LOG(info) << "CbmFsdGeoHandler has initialized maps";
LOG(info) << "fUnitIdToSpecsMap.size() = " << fUnitIdToSpecsMap.size();
LOG(info) << "fModuleIdToSpecsMap.size() = " << fModuleIdToSpecsMap.size();
}
CbmFsdModuleSpecs* CbmFsdGeoHandler::GetModuleSpecsById(Int_t modId)
{
std::map<Int_t, CbmFsdModuleSpecs*>::iterator it;
it = fModuleIdToSpecsMap.find(modId);
if (it == fModuleIdToSpecsMap.end()) return nullptr;
return it->second;
}
CbmFsdUnitSpecs* CbmFsdGeoHandler::GetUnitSpecsById(Int_t unitId)
{
std::map<Int_t, CbmFsdUnitSpecs*>::iterator it;
it = fUnitIdToSpecsMap.find(unitId);
if (it == fUnitIdToSpecsMap.end()) return nullptr;
return it->second;
}
int32_t CbmFsdGeoHandler::GetAddress(TString geoPath) const
{
Int_t moduleID = GetCopyNumberByKey(geoPath, fModuleStr);
Int_t unitID = GetCopyNumberByKey(geoPath, fUnitStr);
return CbmFsdAddress::GetAddress(unitID, moduleID);
}
int32_t CbmFsdGeoHandler::GetCurrentAddress(TVirtualMC* vmc) const
{
Int_t moduleID = -1;
Int_t unitID = -1;
Int_t upstreamOffset = 0;
while (((TString) vmc->CurrentVolOffName(upstreamOffset)).Length() > 0) {
if (((TString) vmc->CurrentVolOffName(upstreamOffset)).Contains(fUnitStr))
vmc->CurrentVolOffID(upstreamOffset, unitID);
if (((TString) vmc->CurrentVolOffName(upstreamOffset)).Contains(fModuleStr))
vmc->CurrentVolOffID(upstreamOffset, moduleID);
upstreamOffset++;
}
return CbmFsdAddress::GetAddress(unitID, moduleID);
}
int32_t CbmFsdGeoHandler::GetCurrentAddress(TGeoManager* geoMan) const
{
Int_t moduleID = -1;
Int_t unitID = -1;
Int_t upstreamOffset = 0;
while (upstreamOffset <= geoMan->GetLevel()) {
TGeoNode* node = geoMan->GetMother(upstreamOffset);
if (((TString) node->GetVolume()->GetName()).Contains(fUnitStr)) unitID = node->GetNumber();
if (((TString) node->GetVolume()->GetName()).Contains(fModuleStr)) moduleID = node->GetNumber();
upstreamOffset++;
}
return CbmFsdAddress::GetAddress(unitID, moduleID);
}
Int_t CbmFsdGeoHandler::GetCopyNumberByKey(TString geoPath, TString key) const
{
Int_t copyNum = -1;
// sanity checks
if (!geoPath.Contains(fBranchStr)) { LOG(warning) << __func__ << ": In geoPath " << fBranchStr << " was not found!"; }
else if (!geoPath.Contains(key)) {
LOG(warning) << __func__ << ": In geoPath " << key << " was not found!";
}
else {
Ssiz_t keyStart = geoPath.Index(key);
TString keyName = geoPath(keyStart, geoPath.Index("/", keyStart) - keyStart);
Ssiz_t copyNumStart = keyName.Last('_') + 1;
TString copyNumStr = keyName(copyNumStart, keyName.Length() - copyNumStart);
if (!copyNumStr.IsDigit()) {
LOG(warning) << __func__ << ": Expected numerical part from " << geoPath << " using key " << key << " is "
<< copyNumStr << " which does not contain only digits!";
}
else {
copyNum = copyNumStr.Atoi();
}
}
return copyNum;
}
/* Copyright (C) 2023 Physikalisches Institut, Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Lukas Chlad [committer] */
#ifndef CBMFSDGEOHANDLER_H
#define CBMFSDGEOHANDLER_H
#include <RtypesCore.h> // for Int_t, Bool_t
#include <TString.h> // for TString
#include <iostream> // for string
#include <map> // for map
#include <vector> // for vector
struct CbmFsdModuleSpecs;
struct CbmFsdUnitSpecs;
class TVirtualMC;
class TGeoManager;
class CbmFsdGeoHandler {
private:
/** Default constructor **/
CbmFsdGeoHandler();
/** Default destructor **/
~CbmFsdGeoHandler() = default;
/** @brief Helper function to extract copy number from geoPath using key word
** @return integer corresponding to copy number of the key word
**/
Int_t GetCopyNumberByKey(TString geoPath, TString key) const;
public:
CbmFsdGeoHandler(const CbmFsdGeoHandler&) = delete;
CbmFsdGeoHandler& operator=(const CbmFsdGeoHandler&) = delete;
/**
* Return Instance of CbmFsdGeoHandler.
*/
static CbmFsdGeoHandler& GetInstance()
{
static CbmFsdGeoHandler fInstance;
return fInstance;
}
/*
* \brief Return CbmFsdModuleSpecs by digi moduleId.
*/
CbmFsdModuleSpecs* GetModuleSpecsById(Int_t id);
/*
* \brief Return CbmFsdUnitSpecs by digi unitId.
*/
CbmFsdUnitSpecs* GetUnitSpecsById(Int_t id);
/*
* \brief Initialize maps.
*/
void InitMaps();
/** @brief Get the unique address from geometry path string
** @return integer corresponding to CbmFsdAddress scheme
**/
int32_t GetAddress(TString geoPath) const;
/** @brief Get the unique address from TVirtualMC
** @return integer corresponding to CbmFsdAddress scheme
**/
int32_t GetCurrentAddress(TVirtualMC* vmc) const;
/** @brief Get the unique address from TGeoManager
** @return integer corresponding to CbmFsdAddress scheme
**/
int32_t GetCurrentAddress(TGeoManager* geoMan) const;
private:
std::map<Int_t, CbmFsdModuleSpecs*> fModuleIdToSpecsMap;
std::map<Int_t, CbmFsdUnitSpecs*> fUnitIdToSpecsMap;
// these key TStrings must corresponds to what one introduce in create geo macros!!
const TString fBranchStr = "fsd_";
const TString fUnitStr = "unit";
const TString fModuleStr = "module";
const TString fActiveMatStr = "scint";
};
#endif /* CBMFSDGEOHANDLER_H */
......@@ -5,6 +5,7 @@ add_subdirectory(bmon)
add_subdirectory(much)
add_subdirectory(mvd)
add_subdirectory(psd)
add_subdirectory(fsd)
add_subdirectory(rich)
add_subdirectory(sts)
add_subdirectory(tof)
......
......
set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}
)
set(SRCS
CbmFsdMC.cxx
CbmFsdDigitize.cxx
)
set(LIBRARY_NAME CbmFsdSim)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
set(PUBLIC_DEPENDENCIES
CbmBase
CbmData
FairRoot::Base
ROOT::Core
ROOT::Physics
)
set(PRIVATE_DEPENDENCIES
CbmSimBase
CbmFsdBase
FairRoot::ParBase
FairLogger::FairLogger
ROOT::Geom
ROOT::MathCore
${VMCLIB}
)
generate_cbm_library()
/* Copyright (C) 2023 Physikalisches Institut Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Alla Maevskaya, Selim Seddiki, Sergey Morozov, Volker Friese, Evgeny Kashirin, Lukas Chlad [committer] */
#include "CbmFsdDigitize.h"
#include "CbmFsdAddress.h"
#include "CbmFsdDigi.h"
#include "CbmFsdDigiPar.h"
#include "CbmFsdPoint.h"
#include "CbmLink.h"
#include "CbmMatch.h"
#include <FairRootManager.h>
#include <FairRun.h>
#include <FairRuntimeDb.h>
#include <Logger.h>
#include <TArrayD.h>
#include <TClonesArray.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TStopwatch.h>
#include <cassert>
#include <iomanip>
#include <iostream>
#include <map>
using std::cout;
using std::endl;
using std::fixed;
using std::right;
using std::setprecision;
using std::setw;
using std::map;
using std::pair;
// ----- Public method Init --------------------------------------------
InitStatus CbmFsdDigitize::Init()
{
if (!fEventMode) { LOG(info) << GetName() << " uses TimeBased mode."; }
else {
LOG(info) << GetName() << " uses Events mode.";
}
// Get RootManager
FairRootManager* ioman = FairRootManager::Instance();
// Get input array
fPointArray = dynamic_cast<TClonesArray*>(ioman->GetObject("FsdPoint"));
if (nullptr == fPointArray) {
LOG(error) << GetName() << ": Error in reading from file! No FsdPoint array.";
return kERROR;
}
TString objectName = fPointArray->GetClass()->GetName();
if (0 != objectName.CompareTo("CbmFsdPoint")) {
LOG(fatal) << GetName() << ": TClonesArray does not contain data of the expected class CbmFsdPoint";
}
// Initialise parameters
InitParams();
// Create and register output array
RegisterOutput();
// Statistics
fNumEvents = 0;
fNumPoints = 0.;
fNumDigis = 0.;
fTimeTot = 0.;
LOG(info) << GetName() << ": Initialisation successful " << kSUCCESS;
return kSUCCESS;
}
// -------------------------------------------------------------------------
void CbmFsdDigitize::SetParContainers()
{
LOG(info) << GetName() << ": Get the digi parameters for FSD";
// Get Base Container
FairRun* ana = FairRun::Instance();
FairRuntimeDb* rtdb = ana->GetRuntimeDb();
fDigiPar = dynamic_cast<CbmFsdDigiPar*>(rtdb->getContainer("CbmFsdDigiPar"));
if (fDigiPar == nullptr) {
LOG(fatal) << GetName() << ": parameter container CbmFsdDigiPar not available in current RuntimeDb!!";
}
}
void CbmFsdDigitize::InitParams()
{
if (!fDigiPar) { LOG(fatal) << GetName() << ": parameter container CbmFsdDigiPar not found!!"; }
// Get values from loaded parameter container
fNumPhotoDets = fDigiPar->GetNumPhotoDets();
fNumUnits = fDigiPar->GetNumUnits();
if (fNumPhotoDets < 1 || fNumUnits < 1) {
LOG(fatal) << GetName() << ": parameter values for FSD digitization does not meet a sanity check!!";
}
fTimeResolution.Set(fNumUnits);
fEnergyResolution.Set(fNumUnits);
fDeadTime.Set(fNumUnits);
for (Int_t iu = 0; iu < fNumUnits; iu++) {
fTimeResolution[iu] = fDigiPar->GetTimeResolution(iu);
fEnergyResolution[iu] = fDigiPar->GetEnergyResolution(iu);
fDeadTime[iu] = fDigiPar->GetDeadTime(iu);
if (fDeadTime[iu] < 0. || fTimeResolution[iu] < 0. || fEnergyResolution[iu] < 0.) {
LOG(fatal) << GetName() << ": parameter values for FSD digitization does not meet a sanity check!!";
}
}
}
// ----- Public method Exec --------------------------------------------
void CbmFsdDigitize::Exec(Option_t*)
{
// Event info (for event time)
GetEventInfo();
TStopwatch timer;
timer.Start();
Int_t nDigis = 0;
// digis that are distant in time to this event (and thus can not be modified anymore) should be send to DAQ
if (!fEventMode) ReleaseBuffer(kFALSE);
Int_t nPoints = fPointArray->GetEntriesFast();
LOG(debug) << fName << ": processing event " << fCurrentEvent << " at t = " << fCurrentEventTime << " ns";
// Declare some variables
CbmFsdPoint* point = nullptr;
Int_t modId = -1;
Int_t unitId = -1;
Int_t photodetId = -1;
if (fNumPhotoDets > 1) {
// random uniform split into photo detectors?!
photodetId = gRandom->Integer(static_cast<UInt_t>(fNumPhotoDets));
}
else {
photodetId = 0;
}
// Loop over FsdPoints
for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
point = static_cast<CbmFsdPoint*>(fPointArray->At(iPoint));
if (!point) continue;
int32_t pointAddress = static_cast<int32_t>(point->GetDetectorID());
modId =
static_cast<Int_t>(CbmFsdAddress::GetElementId(pointAddress, static_cast<int32_t>(CbmFsdAddress::Level::Module)));
unitId =
static_cast<Int_t>(CbmFsdAddress::GetElementId(pointAddress, static_cast<int32_t>(CbmFsdAddress::Level::Unit)));
Double_t eloss = point->GetEnergyLoss() + gRandom->Gaus(0., fEnergyResolution[unitId]);
Double_t time = fCurrentEventTime + point->GetTime() + gRandom->Gaus(0., fTimeResolution[unitId]);
uint32_t address = static_cast<uint32_t>(CbmFsdAddress::GetAddress(
static_cast<uint32_t>(unitId), static_cast<uint32_t>(modId), static_cast<uint32_t>(photodetId)));
CbmLink link(eloss, iPoint, fCurrentMCEntry, fCurrentInput); // weight of link is the energy loss
auto it = fDigiBuffer.find(pointAddress);
if (it != fDigiBuffer.end()) {
// there is already in buffer a digi with this key
Double_t timeDiff = std::fabs(time - it->second.first->GetTime());
if (timeDiff < fDeadTime[unitId]) {
// modify found digi : add energy deposition, modify time if earlier
it->second.first->SetEdep(it->second.first->GetEdep() + eloss);
if (time < it->second.first->GetTime()) it->second.first->SetTime(time);
// add link to digimatch
it->second.second->AddLink(link);
}
else {
// not within time cut -> send the digi from buffer and replace by the new one
CbmFsdDigi* oldDigi =
new CbmFsdDigi(it->second.first->GetAddress(), it->second.first->GetTime(), it->second.first->GetEdep());
CbmMatch* oldDigiMatch = new CbmMatch(*it->second.second);
if (fCreateMatches) SendData(oldDigi->GetTime(), oldDigi, oldDigiMatch);
else
SendData(oldDigi->GetTime(), oldDigi);
fDigiBuffer.erase(it);
CbmFsdDigi* digi = new CbmFsdDigi(address, time, eloss);
CbmMatch* match = new CbmMatch();
match->AddLink(link);
fDigiBuffer.insert(std::make_pair(pointAddress, std::make_pair(digi, match)));
nDigis++;
}
}
else {
// digi with this key is not yet in buffer -> insert it
CbmFsdDigi* digi = new CbmFsdDigi(address, time, eloss);
CbmMatch* match = new CbmMatch();
match->AddLink(link);
fDigiBuffer.insert(std::make_pair(pointAddress, std::make_pair(digi, match)));
nDigis++;
}
} // Loop over MCPoints
// in EventMode send and clear the whole buffer after MCEvent is processed
if (fEventMode) ReleaseBuffer(kTRUE);
// --- Event log
timer.Stop();
LOG(info) << "+ " << setw(15) << GetName() << ": Event " << setw(6) << right << fCurrentEvent << " at " << fixed
<< setprecision(3) << fCurrentEventTime << " ns, points: " << nPoints << ", digis: " << nDigis
<< ". Exec time " << setprecision(6) << timer.RealTime() << " s.";
// --- Run statistics
fNumEvents++;
fNumPoints += nPoints;
fNumDigis += nDigis;
fTimeTot += timer.RealTime();
}
// -------------------------------------------------------------------------
// ----- End-of-run ----------------------------------------------------
void CbmFsdDigitize::Finish()
{
std::cout << std::endl;
LOG(info) << "=====================================";
LOG(info) << GetName() << ": Finish run";
if (!fEventMode) ReleaseBuffer(kTRUE); // in time-based mode, send all what is left in buffer to DAQ
LOG(info) << GetName() << ": Run summary";
LOG(info) << "Events processed : " << fNumEvents;
LOG(info) << "FsdPoint / event : " << setprecision(1) << fNumPoints / Double_t(fNumEvents);
LOG(info) << "FsdDigi / event : " << fNumDigis / Double_t(fNumEvents);
LOG(info) << "Digis per point : " << setprecision(6) << fNumDigis / fNumPoints;
LOG(info) << "Real time per event : " << fTimeTot / Double_t(fNumEvents) << " s";
LOG(info) << "=====================================";
}
// -------------------------------------------------------------------------
void CbmFsdDigitize::ReleaseBuffer(Bool_t sendEverything)
{
if (sendEverything) {
// send everything
for (const auto& dp : fDigiBuffer) {
CbmFsdDigi* digi =
new CbmFsdDigi(dp.second.first->GetAddress(), dp.second.first->GetTime(), dp.second.first->GetEdep());
CbmMatch* digiMatch = new CbmMatch(*dp.second.second);
if (fCreateMatches) SendData(digi->GetTime(), digi, digiMatch);
else
SendData(digi->GetTime(), digi);
} // # digi buffer
fDigiBuffer.clear();
}
else {
// send only if time difference to the new start of event is larger than deadtime
std::vector<int32_t> keysSent;
for (const auto& dp : fDigiBuffer) {
Int_t unitId = static_cast<Int_t>(CbmFsdAddress::GetElementId(static_cast<int32_t>(dp.second.first->GetAddress()),
static_cast<int32_t>(CbmFsdAddress::Level::Unit)));
Double_t timeDiff = std::fabs(dp.second.first->GetTime() - fCurrentEventTime);
if (timeDiff > fDeadTime[unitId]) {
// send this digi
CbmFsdDigi* digi =
new CbmFsdDigi(dp.second.first->GetAddress(), dp.second.first->GetTime(), dp.second.first->GetEdep());
CbmMatch* digiMatch = new CbmMatch(*dp.second.second);
if (fCreateMatches) SendData(digi->GetTime(), digi, digiMatch);
else
SendData(digi->GetTime(), digi);
// save which keys were sent and should be removed
keysSent.push_back(dp.first);
} // ? time
} // # digi buffer
// remove the sent elements from buffer
for (auto& key : keysSent) {
fDigiBuffer.erase(key);
}
} // ? sendEverything
}
ClassImp(CbmFsdDigitize)
/* Copyright (C) 2023 Physikalisches Institut Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Sergey Morozov, Volker Friese, Lukas Chlad [committer] */
/** @class CbmFsdDigitize
** @date 14.07.2023
** @author Lukas Chlad <l.chlad@gsi.de>
** @brief Class for the digitization of the CBM-FSD
**
** The digitizer produces digits of type CbmFsdDigi as sum of Edep of Points and adds smearing in energy and time according to parameters
**/
#ifndef CBMFSDDIGITIZE_H
#define CBMFSDDIGITIZE_H 1
#include "CbmDefs.h"
#include "CbmDigitize.h"
#include "CbmFsdDigi.h"
#include <TArrayD.h>
class TClonesArray;
class CbmFsdDigiPar;
class CbmFsdDigitize : public CbmDigitize<CbmFsdDigi> {
public:
/** Default constructor **/
CbmFsdDigitize() : CbmDigitize<CbmFsdDigi>("FsdDigitize") {};
/** Destructor **/
virtual ~CbmFsdDigitize() = default;
CbmFsdDigitize(const CbmFsdDigitize&) = delete;
CbmFsdDigitize operator=(const CbmFsdDigitize&) = delete;
ECbmModuleId GetSystemId() const { return ECbmModuleId::kFsd; }
/**
** @brief Inherited from FairTask.
**/
virtual InitStatus Init();
/**
** @brief Inherited from FairTask.
**/
virtual void SetParContainers();
/** Virtual method Exec **/
virtual void Exec(Option_t* opt);
/** @brief End-of-run action **/
virtual void Finish();
private:
CbmFsdDigiPar* fDigiPar = nullptr;
Int_t fNumPhotoDets = -1;
Int_t fNumUnits = -1;
TArrayD fTimeResolution {};
TArrayD fEnergyResolution {};
TArrayD fDeadTime {};
Int_t fNumEvents = 0;
Double_t fNumPoints = 0.;
Double_t fNumDigis = 0.;
Double_t fTimeTot = 0.;
/** Input array of CbmFsdPoints **/
TClonesArray* fPointArray = nullptr;
// Temporary storage for digis, key is DetectorID from FsdPoint
std::map<int32_t, std::pair<CbmFsdDigi*, CbmMatch*>> fDigiBuffer;
/** @brief Initialise the parameters **/
void InitParams();
/** @brief release digi from local buffer to CbmDaq
** TIME BASED
** - at the beginning of MCEvent loop over digi buffer
** and send those digis that are too far from start of current event to be possibly edited
** free the location of sent digis
** - at the end of whole run do the same only send everything, clear buffer
**
** EVENT BASED
** - at the beginning of MCEvent loop send everything, clear buffer
**/
void ReleaseBuffer(Bool_t sendEverything);
ClassDef(CbmFsdDigitize, 1);
};
#endif
/* Copyright (C) 2023 Physikalisches Institut Eberhard Karls Universitaet Tuebingen, Tuebingen
SPDX-License-Identifier: GPL-3.0-only
Authors: Alla Maevskaya, Florian Uhlig, Lukas Chlad [committer] */
/** @file CbmFsdMC.cxx
** @since 14.07.2023
** @author Lukas Chlad <l.chlad@gsi.de>
**
**/
#include "CbmFsdMC.h"
#include "CbmFsdGeoHandler.h"
#include "CbmFsdPoint.h"
#include "CbmGeometryUtils.h"
#include "CbmModuleList.h"
#include "CbmStack.h"
#include <FairVolume.h>
#include <TGeoNode.h>
#include <TGeoVolume.h>
#include <TVirtualMC.h>
#include <cassert>
#include <string>
// ----- Destructor ----------------------------------------------------
CbmFsdMC::~CbmFsdMC()
{
if (fFsdPoints) {
fFsdPoints->Delete();
delete fFsdPoints;
}
}
// -------------------------------------------------------------------------
// ----- Construct the geometry from file ------------------------------
void CbmFsdMC::ConstructGeometry()
{
LOG(info) << "Importing FSD geometry from ROOT file " << fgeoName.Data();
Cbm::GeometryUtils::ImportRootGeometry(fgeoName, this);
}
// -------------------------------------------------------------------------
// ----- End of event action -------------------------------------------
void CbmFsdMC::EndOfEvent()
{
Print(); // Status output
fFsdPoints->Delete();
}
// -------------------------------------------------------------------------
// ----- Print ---------------------------------------------------------
void CbmFsdMC::Print(Option_t*) const
{
LOG(info) << fName << ": " << fFsdPoints->GetEntriesFast() << " points registered in this event.";
}
// -------------------------------------------------------------------------
// ----- Initialise ----------------------------------------------------
void CbmFsdMC::Initialize()
{
// --- Instantiate the output array
fFsdPoints = new TClonesArray("CbmFsdPoint");
// --- Call the Initialise method of the mother class
FairDetector::Initialize();
}
// -------------------------------------------------------------------------
// ----- Public method ProcessHits --------------------------------------
Bool_t CbmFsdMC::ProcessHits(FairVolume*)
{
// No action for neutral particles
if (TMath::Abs(gMC->TrackCharge()) <= 0) return kFALSE;
// --- If this is the first step for the track in the volume:
// Reset energy loss and store track parameters
if (gMC->IsTrackEntering()) {
fTrackID = gMC->GetStack()->GetCurrentTrackNumber();
fAddress = CbmFsdGeoHandler::GetInstance().GetCurrentAddress(gMC);
gMC->TrackPosition(fPos);
gMC->TrackMomentum(fMom);
fTime = gMC->TrackTime() * 1.0e09;
fLength = gMC->TrackLength();
fEloss = 0.;
} //? track entering
// --- For all steps within active volume: sum up differential energy loss
fEloss += gMC->Edep();
// --- If track is leaving: get track parameters and create CbmstsPoint
if (gMC->IsTrackExiting() || gMC->IsTrackStop() || gMC->IsTrackDisappeared()) {
// Create CbmFsdPoint
Int_t size = fFsdPoints->GetEntriesFast();
new ((*fFsdPoints)[size]) CbmFsdPoint(fTrackID, fAddress, fPos.Vect(), fMom.Vect(), fTime, fLength, fEloss);
// --- Increment number of FsdPoints for this track in the stack
CbmStack* stack = dynamic_cast<CbmStack*>(gMC->GetStack());
assert(stack);
stack->AddPoint(ECbmModuleId::kFsd);
} //? track is exiting or stopped
return kTRUE;
}
// -------------------------------------------------------------------------
// ----- Register the sensitive volumes --------------------------------
void CbmFsdMC::RegisterSensitiveVolumes(TGeoNode* node)
{
TObjArray* daughters = node->GetVolume()->GetNodes();
for (Int_t iDaughter = 0; iDaughter < daughters->GetEntriesFast(); iDaughter++) {
TGeoNode* daughter = dynamic_cast<TGeoNode*>(daughters->At(iDaughter));
assert(daughter);
if (daughter->GetNdaughters() > 0) RegisterSensitiveVolumes(daughter);
TGeoVolume* daughterVolume = daughter->GetVolume();
if (CheckIfSensitive(daughterVolume->GetName())) { AddSensitiveVolume(daughterVolume); } //? Sensitive volume
} //# Daughter nodes
}
// -------------------------------------------------------------------------
ClassImp(CbmFsdMC)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment