diff --git a/CMakeLists.txt b/CMakeLists.txt index b54c3eed..b75f8ba0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -267,6 +267,7 @@ add_subdirectory (field) add_subdirectory (genfit) add_subdirectory (millepede) add_subdirectory (analysis) +add_subdirectory (digitization) FIND_PATH(TEvePath NAMES TEveEventManager.h PATHS ${SIMPATH}/tools/root/include diff --git a/digitization/CMakeLists.txt b/digitization/CMakeLists.txt new file mode 100644 index 00000000..4a7fd9a2 --- /dev/null +++ b/digitization/CMakeLists.txt @@ -0,0 +1,31 @@ +set(INCLUDE_DIRECTORIES +${BASE_INCLUDE_DIRECTORIES} +${CMAKE_SOURCE_DIR}/shipdata +${CMAKE_SOURCE_DIR}/shipLHC +${CMAKE_SOURCE_DIR}/digitization +) + +include_directories( ${INCLUDE_DIRECTORIES} ${VMC_INCLUDE_DIRS} ${FAIRROOT_INCLUDE_DIR} ${FairLogger_INCDIR} ${FMT_INCLUDE_DIR} +) +include_directories(SYSTEM ${SYSTEM_INCLUDE_DIRECTORIES}) + +set(LINK_DIRECTORIES +${ROOT_LIBRARY_DIR} +${FAIRROOT_LIBRARY_DIR} +) + +link_directories( ${LINK_DIRECTORIES}) + +set(SRCS +digitize.cxx +) + +Set(HEADERS) +Set(LINKDEF digitizationLinkDef.h) +Set(LIBRARY_NAME digitization) +Set(DEPENDENCIES Base ShipData shipLHC GeoBase ParBase Geom Core) + +GENERATE_LIBRARY() + +add_executable(digitize_mc digitize_mc.cxx) +target_link_libraries(digitize_mc ${ROOT_LIBRARIES} digitization veto FairTools -lROOTTPython) diff --git a/digitization/digitizationLinkDef.h b/digitization/digitizationLinkDef.h new file mode 100644 index 00000000..2d33e52d --- /dev/null +++ b/digitization/digitizationLinkDef.h @@ -0,0 +1,24 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ nestedclasses; +#pragma link C++ nestedtypedef; +#pragma link C++ namespace advsnd; +#pragma link C++ defined_in namespace advsnd; + +#pragma link C++ class advsnd::DigitizePoints+; +#pragma link C++ class advsnd::DigitizePoints+; +#pragma link C++ class advsnd::LinkPointsToDigi+; +#pragma link C++ class advsnd::LinkPointsToDigi+; + +#pragma link C++ function advsnd::GetStripId; +#pragma link C++ function advsnd::GetStripId; +#pragma link C++ function advsnd::Digitize; +#pragma link C++ function advsnd::Digitize; +#pragma link C++ function advsnd::McLink; +#pragma link C++ function advsnd::McLink; + +#endif \ No newline at end of file diff --git a/digitization/digitize.cxx b/digitization/digitize.cxx new file mode 100644 index 00000000..90d445d7 --- /dev/null +++ b/digitization/digitize.cxx @@ -0,0 +1 @@ +#include "digitize.h" \ No newline at end of file diff --git a/digitization/digitize.h b/digitization/digitize.h new file mode 100644 index 00000000..70496fe9 --- /dev/null +++ b/digitization/digitize.h @@ -0,0 +1,118 @@ +#ifndef ADVSND_DIGITIZE_H +#define ADVSND_DIGITIZE_H + +#include "TGeoNavigator.h" +#include "AdvTargetPoint.h" +#include "AdvTargetHit.h" +#include "AdvMuFilterPoint.h" +#include "AdvMuFilterHit.h" +#include "Hit2MCPoints.h" +#include "TString.h" + +#include +#include +#include +#include +#include + +namespace advsnd { + + template + int GetStripId(const T& point, TGeoNavigator* nav) { + auto detID = point.GetDetectorID(); + TString path; + + if constexpr (std::is_same::value) { + path = TString::Format( + "/cave_1/Detector_0/volAdvTarget_0/Target_Layer_%d/SensorModule_%d/Target_SensorVolume_%d", + point.GetLayer(), point.GetModule(), detID); + } + else if constexpr (std::is_same::value) { + path = TString::Format( + "/cave_1/Detector_0/volAdvMuFilter_0/HCAL_Layer_%d/SensorModule_%d/HCAL_SensorVolume_%d", + point.GetLayer(), point.GetModule(), detID); + } + else { + std::cerr << "Invalid type\n"; + } + + if (nav->CheckPath(path)) { + nav->cd(path); + } else { + std::cerr << "Invalid geometry path: " << path << "\n"; + } + + nav->GetCurrentNode(); + double global_pos[3] = {point.GetX(), point.GetY(), point.GetZ()}; + double local_pos[3]; + nav->MasterToLocal(global_pos, local_pos); + + int strip = floor((local_pos[1] / (advsnd::sensor_length / advsnd::strips)) + (advsnd::strips / 2)); + strip = std::clamp(strip, 0, advsnd::strips - 1); + return strip; + } + + template + std::vector Digitize(const std::vector& v_points, TGeoNavigator* nav) { + std::map> hit_collector{}; + std::vector v_hits; + + // --- Loop over all target points --- + for (const auto& point : v_points) { + int strip = GetStripId(point, nav); + auto detector_id = point.GetDetectorID() - 999 + strip; + hit_collector[detector_id].emplace_back(&point); + } + // --- Build hits --- + for (const auto& [detector_id, points] : hit_collector) { + v_hits.emplace_back(U(detector_id, points)); + } + return v_hits; + } + + template + Hit2MCPoints McLink(const std::vector& v_points, TGeoNavigator* nav) { + int point_index{0}; + std::vector v_detector_id; + Hit2MCPoints mc_links; + std::map> mc_points{}; + std::map norm{}; + + // --- Loop over all points --- + for (const auto& point : v_points) { + int strip = GetStripId(point, nav); + const int detector_id = point.GetDetectorID() - 999 + strip; + v_detector_id.emplace_back(detector_id); + mc_points[detector_id][point_index++] = point.GetEnergyLoss(); + norm[detector_id] += point.GetEnergyLoss(); + } + + // --- Build MC links --- + for (const auto detector_id : v_detector_id) { + for (const auto& [point_id, energy_loss] : mc_points[detector_id]) { + mc_links.Add(detector_id, point_id, energy_loss / norm[detector_id]); + } + } + return mc_links; + } + + template + class DigitizePoints { + public: + DigitizePoints(TGeoNavigator* nav) : m_nav(nav) {} + std::vector operator()(const std::vector& v_points) const { return Digitize(v_points, m_nav); } + private: + TGeoNavigator* m_nav; + }; + + template + class LinkPointsToDigi { + public: + LinkPointsToDigi(TGeoNavigator* nav) : m_nav(nav) {} + Hit2MCPoints operator()(const std::vector& v_points) const { return McLink(v_points, m_nav); } + private: + TGeoNavigator* m_nav; + }; +} + +#endif \ No newline at end of file diff --git a/digitization/digitize_mc.cxx b/digitization/digitize_mc.cxx new file mode 100644 index 00000000..d526358c --- /dev/null +++ b/digitization/digitize_mc.cxx @@ -0,0 +1,45 @@ +#include "TPython.h" +#include "TGeoNavigator.h" +#include "TGeoManager.h" +#include "AdvTargetPoint.h" +#include "AdvTargetHit.h" +#include "AdvMuFilterPoint.h" +#include "AdvMuFilterHit.h" +#include "Hit2MCPoints.h" +//#include "FairLogger.h" +#include "ROOT/RDataFrame.hxx" +#include "digitize.h" + +#include +#include + +TGeoNavigator* initGeometry(const std::string& geometry_path) +{ + TPython::Exec("import SndlhcGeo"); + TPython::Exec(("SndlhcGeo.GeoInterface('" + geometry_path + "')").c_str()); + + if (!gGeoManager) { + std::cerr << "Geofile required\n"; + } + TGeoNavigator* nav = gGeoManager->GetCurrentNavigator(); + if (!nav) { + std::cerr << "Failed to get TGeoNavigator from gGeoManager.\n"; + } + return nav; +} + +int main() { + ROOT::RDF::RSnapshotOptions opts; + opts.fOutputFormat = ROOT::RDF::ESnapshotOutputFormat::kRNTuple; + + auto df = ROOT::RDataFrame("cbmsim", "/afs/cern.ch/work/f/fmei/private/RNTuples-for-advsndsw/data/rntuples-sndLHC.Ntuple-TGeant4.root"); + + std::string geometry_path = "/afs/cern.ch/work/f/fmei/private/RNTuples-for-advsndsw/data_multi/geofile_full.Ntuple-TGeant4.root"; + TGeoNavigator* nav = initGeometry(geometry_path); + + auto df2 = df.Define("Digi_AdvTargetHits", advsnd::DigitizePoints(nav), {"AdvTargetPoint"}); + auto df3 = df2.Define("Digi_AdvTargetHits2MCPoints", advsnd::LinkPointsToDigi(nav), {"AdvTargetPoint"}); + auto df4 = df3.Define("Digi_AdvMuFilterHits", advsnd::DigitizePoints(nav), {"AdvMuFilterPoint"}); + auto df5 = df4.Define("Digi_AdvMuFilterHits2MCPoints", advsnd::LinkPointsToDigi(nav), {"AdvMuFilterPoint"}); + df5.Snapshot("cbmsim", "/afs/cern.ch/work/f/fmei/private/RNTuples-for-advsndsw/data_multi/rntuples-df-compiled-sndLHC.Ntuple-TGeant4_digCPP.root", {"AdvMuFilterPoint", "AdvTargetPoint", "Digi_AdvTargetHits", "Digi_AdvTargetHits2MCPoints", "Digi_AdvMuFilterHits","Digi_AdvMuFilterHits2MCPoints"}, opts); +} \ No newline at end of file diff --git a/shipLHC/AdvMuFilterHit.cxx b/shipLHC/AdvMuFilterHit.cxx index 23b570a5..d436abeb 100644 --- a/shipLHC/AdvMuFilterHit.cxx +++ b/shipLHC/AdvMuFilterHit.cxx @@ -36,6 +36,16 @@ AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector LOG(DEBUG) << "signal created"; } +AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector& V) + : SndlhcHit(detID) +{ + flag = true; + for (Int_t i = 0; i < 16; i++) { + fMasked[i] = kFALSE; + } + LOG(DEBUG) << "signal created"; +} + // ----- Public method Print ------------------------------------------- void AdvMuFilterHit::Print() const { diff --git a/shipLHC/AdvMuFilterHit.h b/shipLHC/AdvMuFilterHit.h index 8497ae29..b9495a75 100644 --- a/shipLHC/AdvMuFilterHit.h +++ b/shipLHC/AdvMuFilterHit.h @@ -15,6 +15,7 @@ class AdvMuFilterHit : public SndlhcHit // Constructor from MuFilterPoint AdvMuFilterHit(Int_t detID, const std::vector&); + AdvMuFilterHit(Int_t detID, const std::vector&); /** Destructor **/ ~AdvMuFilterHit() = default; @@ -24,14 +25,14 @@ class AdvMuFilterHit : public SndlhcHit bool isValid() const { return flag; } bool isMasked(Int_t i) const { return fMasked[i]; } void SetMasked(Int_t i) { fMasked[i] = kTRUE; } - int constexpr GetLayer() { return fDetectorID >> 17; } - int constexpr GetPlane() { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane - int constexpr GetRow() { return (fDetectorID >> 13) % 8; } - int constexpr GetColumn() { return (fDetectorID >> 11) % 4; } - int constexpr GetSensor() { return (fDetectorID >> 10) % 2; } - int constexpr GetStrip() { return (fDetectorID) % 1024; } - int constexpr GetModule() { return advsnd::hcal::columns * GetRow() + 1 + GetColumn(); } - bool constexpr isVertical() { return GetPlane() == 1; }; + int constexpr GetLayer() const { return fDetectorID >> 17; } + int constexpr GetPlane() const { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane + int constexpr GetRow() const { return (fDetectorID >> 13) % 8; } + int constexpr GetColumn() const { return (fDetectorID >> 11) % 4; } + int constexpr GetSensor() const { return (fDetectorID >> 10) % 2; } + int constexpr GetStrip() const { return (fDetectorID) % 1024; } + int constexpr GetModule() const { return advsnd::hcal::columns * GetRow() + 1 + GetColumn(); } + bool constexpr isVertical() const { return GetPlane() == 1; }; private: bool flag; ///< flag diff --git a/shipLHC/AdvMuFilterPoint.h b/shipLHC/AdvMuFilterPoint.h index 57cc09fb..6a212f2f 100644 --- a/shipLHC/AdvMuFilterPoint.h +++ b/shipLHC/AdvMuFilterPoint.h @@ -39,21 +39,18 @@ class AdvMuFilterPoint : public FairMCPoint virtual void Print(const Option_t* opt) const; Int_t PdgCode() const { return fPdgCode; } - int constexpr GetLayer() { return fDetectorID >> 17; } - int constexpr GetPlane() { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane - int constexpr GetRow() { return (fDetectorID >> 13) % 8; } - int constexpr GetColumn() { return (fDetectorID >> 11) % 4; } - int constexpr GetSensor() { return (fDetectorID >> 10) % 2; } - int constexpr GetStrip() { return (fDetectorID) % 1024; } - int constexpr GetModule() { return advsnd::hcal::columns * GetRow() + 1 + GetColumn(); } - bool constexpr isVertical() { return GetPlane() == 1; }; + int constexpr GetLayer() const { return fDetectorID >> 17; } + int constexpr GetPlane() const { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane + int constexpr GetRow() const { return (fDetectorID >> 13) % 8; } + int constexpr GetColumn() const { return (fDetectorID >> 11) % 4; } + int constexpr GetSensor() const { return (fDetectorID >> 10) % 2; } + int constexpr GetStrip() const { return (fDetectorID) % 1024; } + int constexpr GetModule() const { return advsnd::hcal::columns * GetRow() + 1 + GetColumn(); } + bool constexpr isVertical() const { return GetPlane() == 1; }; private: Int_t fPdgCode; - /** Copy constructor **/ - AdvMuFilterPoint(const AdvMuFilterPoint& point); - AdvMuFilterPoint operator=(const AdvMuFilterPoint& point); ClassDef(AdvMuFilterPoint, 1) }; diff --git a/shipLHC/AdvTargetHit.cxx b/shipLHC/AdvTargetHit.cxx index 15686ab5..979ecb6b 100644 --- a/shipLHC/AdvTargetHit.cxx +++ b/shipLHC/AdvTargetHit.cxx @@ -36,6 +36,16 @@ AdvTargetHit::AdvTargetHit(Int_t detID, const std::vector& V) LOG(DEBUG) << "signal created"; } +AdvTargetHit::AdvTargetHit(Int_t detID, const std::vector& V) + : SndlhcHit(detID) +{ + flag = true; + for (Int_t i = 0; i < 16; i++) { + fMasked[i] = kFALSE; + } + LOG(DEBUG) << "signal created"; +} + // ----- Public method Print ------------------------------------------- void AdvTargetHit::Print() const { diff --git a/shipLHC/AdvTargetHit.h b/shipLHC/AdvTargetHit.h index ff8a4917..5f3e92e8 100644 --- a/shipLHC/AdvTargetHit.h +++ b/shipLHC/AdvTargetHit.h @@ -15,6 +15,7 @@ class AdvTargetHit : public SndlhcHit // Constructor from AdvTargetPoint AdvTargetHit(Int_t detID, const std::vector&); + AdvTargetHit(Int_t detID, const std::vector&); /** Destructor **/ ~AdvTargetHit() = default; @@ -24,14 +25,14 @@ class AdvTargetHit : public SndlhcHit bool isValid() const { return flag; } bool isMasked(Int_t i) const { return fMasked[i]; } void SetMasked(Int_t i) { fMasked[i] = kTRUE; } - int constexpr GetLayer() { return fDetectorID >> 17; } - int constexpr GetPlane() { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane - int constexpr GetRow() { return (fDetectorID >> 13) % 8; } - int constexpr GetColumn() { return (fDetectorID >> 11) % 4; } - int constexpr GetSensor() { return (fDetectorID >> 10) % 2; } - int constexpr GetStrip() { return (fDetectorID) % 1024; } - int constexpr GetModule() { return advsnd::target::columns * GetRow() + 1 + GetColumn(); } - bool constexpr isVertical() { return GetPlane() == 1; }; + int constexpr GetLayer() const { return fDetectorID >> 17; } + int constexpr GetPlane() const { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane + int constexpr GetRow() const { return (fDetectorID >> 13) % 8; } + int constexpr GetColumn() const { return (fDetectorID >> 11) % 4; } + int constexpr GetSensor() const { return (fDetectorID >> 10) % 2; } + int constexpr GetStrip() const { return (fDetectorID) % 1024; } + int constexpr GetModule() const { return advsnd::target::columns * GetRow() + 1 + GetColumn(); } + bool constexpr isVertical() const { return GetPlane() == 1; }; private: bool flag; ///< flag diff --git a/shipLHC/AdvTargetPoint.h b/shipLHC/AdvTargetPoint.h index 854d0899..94b819bb 100644 --- a/shipLHC/AdvTargetPoint.h +++ b/shipLHC/AdvTargetPoint.h @@ -42,14 +42,14 @@ class AdvTargetPoint : public FairMCPoint virtual void Print(const Option_t* opt) const; Int_t PdgCode() const { return fPdgCode; } - int constexpr GetLayer() { return fDetectorID >> 17; } - int constexpr GetPlane() { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane - int constexpr GetRow() { return (fDetectorID >> 13) % 8; } - int constexpr GetColumn() { return (fDetectorID >> 11) % 4; } - int constexpr GetSensor() { return (fDetectorID >> 10) % 2; } - int constexpr GetStrip() { return (fDetectorID) % 1024; } - int constexpr GetModule() { return advsnd::target::columns * GetRow() + 1 + GetColumn(); } - bool constexpr IsVertical() { return GetPlane() == 1; }; + int constexpr GetLayer() const { return fDetectorID >> 17; } + int constexpr GetPlane() const { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane + int constexpr GetRow() const { return (fDetectorID >> 13) % 8; } + int constexpr GetColumn() const { return (fDetectorID >> 11) % 4; } + int constexpr GetSensor() const { return (fDetectorID >> 10) % 2; } + int constexpr GetStrip() const { return (fDetectorID) % 1024; } + int constexpr GetModule() const { return advsnd::target::columns * GetRow() + 1 + GetColumn(); } + bool constexpr IsVertical() const { return GetPlane() == 1; }; TVector3 GetEntryPoint() const { return TVector3(2 * fX - fExitX, 2 * fY - fExitY, 2 * fZ - fExitZ); } TVector3 GetExitPoint() const { return TVector3(fExitX, fExitY, fExitZ); } @@ -59,9 +59,6 @@ class AdvTargetPoint : public FairMCPoint Double_t fExitY; Double_t fExitZ; - /** Copy constructor **/ - AdvTargetPoint(const AdvTargetPoint& point); - AdvTargetPoint operator=(const AdvTargetPoint& point); ClassDef(AdvTargetPoint, 1) }; diff --git a/shipdata/SndlhcHit.cxx b/shipdata/SndlhcHit.cxx index f6b6b38e..5843dcc5 100644 --- a/shipdata/SndlhcHit.cxx +++ b/shipdata/SndlhcHit.cxx @@ -25,11 +25,11 @@ SndlhcHit::SndlhcHit(Int_t detID,Int_t nP,Int_t nS) } } -Float_t SndlhcHit::GetSignal(Int_t nChannel) +Float_t SndlhcHit::GetSignal(Int_t nChannel) const { return signals[nChannel]; } -Float_t SndlhcHit::GetTime(Int_t nChannel) +Float_t SndlhcHit::GetTime(Int_t nChannel) const { return times[nChannel]; } diff --git a/shipdata/SndlhcHit.h b/shipdata/SndlhcHit.h index f56ba5ef..672deb74 100644 --- a/shipdata/SndlhcHit.h +++ b/shipdata/SndlhcHit.h @@ -34,18 +34,18 @@ class SndlhcHit : public TObject /** Accessors **/ Int_t GetDetectorID() const { return fDetectorID; }; - Float_t GetSignal(Int_t nChannel=0); - Float_t GetTime(Int_t nChannel=0); + Float_t GetSignal(Int_t nChannel=0) const; + Float_t GetTime(Int_t nChannel=0) const; Int_t GetnSiPMs() const { return nSiPMs; }; Int_t GetnSides() const { return nSides; }; /** Modifiers **/ void SetDigi(Float_t s, Float_t t,Int_t i=0) { signals[i]=trunc(100*s)/100;times[i]=trunc(1000*t)/1000.; } void SetDetectorID(Int_t detID) { fDetectorID = detID; } void SetDaqID(Int_t i, Int_t k, Int_t board_id, Int_t tofpet_id, Int_t tofpet_channel) { fDaqID[i] = k*100000 + board_id * 1000 + tofpet_id * 100 + tofpet_channel; } - Int_t GetBoardID(Int_t i) { return int((fDaqID[i]%100000)/1000);} - Int_t GetTofpetID(Int_t i) { return int((fDaqID[i]%1000)/100);} - Int_t Getchannel(Int_t i) { return fDaqID[i]%100;} - Int_t GetRawHitIndex(Int_t i=0) { return int(fDaqID[i]/100000);} + Int_t GetBoardID(Int_t i) const { return int((fDaqID[i]%100000)/1000);} + Int_t GetTofpetID(Int_t i) const { return int((fDaqID[i]%1000)/100);} + Int_t Getchannel(Int_t i) const { return fDaqID[i]%100;} + Int_t GetRawHitIndex(Int_t i=0) const { return int(fDaqID[i]/100000);} // to be implemented by the subdetector diff --git a/veto/vetoPoint.h b/veto/vetoPoint.h index ace05256..c99d6b5a 100644 --- a/veto/vetoPoint.h +++ b/veto/vetoPoint.h @@ -40,11 +40,8 @@ class vetoPoint : public FairMCPoint TVector3 LastMom() const {return fLmom;} private: - /** Copy constructor **/ - Int_t fPdgCode; TVector3 fLpos,fLmom; - vetoPoint(const vetoPoint& point); - vetoPoint operator=(const vetoPoint& point); + Int_t fPdgCode; ClassDef(vetoPoint,3)