Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
31 changes: 31 additions & 0 deletions digitization/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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)
24 changes: 24 additions & 0 deletions digitization/digitizationLinkDef.h
Original file line number Diff line number Diff line change
@@ -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<AdvTargetPoint, AdvTargetHit>+;
#pragma link C++ class advsnd::DigitizePoints<AdvMuFilterPoint, AdvMuFilterHit>+;
#pragma link C++ class advsnd::LinkPointsToDigi<AdvTargetPoint>+;
#pragma link C++ class advsnd::LinkPointsToDigi<AdvMuFilterPoint>+;

#pragma link C++ function advsnd::GetStripId<AdvTargetPoint>;
#pragma link C++ function advsnd::GetStripId<AdvMuFilterPoint>;
#pragma link C++ function advsnd::Digitize<AdvTargetPoint, AdvTargetHit>;
#pragma link C++ function advsnd::Digitize<AdvMuFilterPoint, AdvMuFilterHit>;
#pragma link C++ function advsnd::McLink<AdvTargetPoint>;
#pragma link C++ function advsnd::McLink<AdvMuFilterPoint>;

#endif
1 change: 1 addition & 0 deletions digitization/digitize.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
#include "digitize.h"
118 changes: 118 additions & 0 deletions digitization/digitize.h
Original file line number Diff line number Diff line change
@@ -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 <vector>
#include <map>
#include <algorithm>
#include <type_traits>
#include <iostream>

namespace advsnd {

template <typename T>
int GetStripId(const T& point, TGeoNavigator* nav) {
auto detID = point.GetDetectorID();
TString path;

if constexpr (std::is_same<T, AdvTargetPoint>::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<T, AdvMuFilterPoint>::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 <typename T, typename U>
std::vector<U> Digitize(const std::vector<T>& v_points, TGeoNavigator* nav) {
std::map<int, std::vector<const T*>> hit_collector{};
std::vector<U> 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 <typename T>
Hit2MCPoints McLink(const std::vector<T>& v_points, TGeoNavigator* nav) {
int point_index{0};
std::vector<int> v_detector_id;
Hit2MCPoints mc_links;
std::map<int, std::map<int, double>> mc_points{};
std::map<int, double> 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 <typename T, typename U>
class DigitizePoints {
public:
DigitizePoints(TGeoNavigator* nav) : m_nav(nav) {}
std::vector<U> operator()(const std::vector<T>& v_points) const { return Digitize<T, U>(v_points, m_nav); }
private:
TGeoNavigator* m_nav;
};

template <typename T>
class LinkPointsToDigi {
public:
LinkPointsToDigi(TGeoNavigator* nav) : m_nav(nav) {}
Hit2MCPoints operator()(const std::vector<T>& v_points) const { return McLink<T>(v_points, m_nav); }
private:
TGeoNavigator* m_nav;
};
}

#endif
45 changes: 45 additions & 0 deletions digitization/digitize_mc.cxx
Original file line number Diff line number Diff line change
@@ -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 <string>
#include <iostream>

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<AdvTargetPoint, AdvTargetHit>(nav), {"AdvTargetPoint"});
auto df3 = df2.Define("Digi_AdvTargetHits2MCPoints", advsnd::LinkPointsToDigi<AdvTargetPoint>(nav), {"AdvTargetPoint"});
auto df4 = df3.Define("Digi_AdvMuFilterHits", advsnd::DigitizePoints<AdvMuFilterPoint, AdvMuFilterHit>(nav), {"AdvMuFilterPoint"});
auto df5 = df4.Define("Digi_AdvMuFilterHits2MCPoints", advsnd::LinkPointsToDigi<AdvMuFilterPoint>(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);
}
10 changes: 10 additions & 0 deletions shipLHC/AdvMuFilterHit.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,16 @@ AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>
LOG(DEBUG) << "signal created";
}

AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector<const AdvMuFilterPoint*>& V)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can safely remove the previous version of this constructor in the lines above
AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>& V)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait a sec, shouldn't this be const std::vector<AdvMuFilterPoint>&?
AFAIK the recommended way to store and read classes in TTree and RNTuple is as std::vector<T> unless not possible.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can safely remove the previous version of this constructor in the lines above AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>& V)

I agree, I did not remove it because i feared it would brake the current digitization

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait a sec, shouldn't this be const std::vector<AdvMuFilterPoint>&? AFAIK the recommended way to store and read classes in TTree and RNTuple is as std::vector<T> unless not possible.

That would be my ideal implementation too, as a minimal change from the current digitization i kept the pointers (also if we implement this as the only constructor the current digitization would definetly break)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we can safely remove the previous version of this constructor in the lines above AdvMuFilterHit::AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>& V)

I agree, I did not remove it because i feared it would brake the current digitization
ok, then we can drop that removal for now, especially since we have a large PR by Nayana that relies on the old format.
We can factorize all these updates.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait a sec, shouldn't this be const std::vector<AdvMuFilterPoint>&? AFAIK the recommended way to store and read classes in TTree and RNTuple is as std::vector<T> unless not possible.

That would be my ideal implementation too, as a minimal change from the current digitization i kept the pointers (also if we implement this as the only constructor the current digitization would definetly break)

hm, then go ahead with your ideal, also advised, implementation, changing also the digitization.
Lets have all that in the PR.

: 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
{
Expand Down
17 changes: 9 additions & 8 deletions shipLHC/AdvMuFilterHit.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ class AdvMuFilterHit : public SndlhcHit

// Constructor from MuFilterPoint
AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>&);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

AdvMuFilterHit(Int_t detID, const std::vector<AdvMuFilterPoint*>&);
this is no longer needed if we have the new version with const pointers

AdvMuFilterHit(Int_t detID, const std::vector<const AdvMuFilterPoint*>&);

/** Destructor **/
~AdvMuFilterHit() = default;
Expand All @@ -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
Expand Down
19 changes: 8 additions & 11 deletions shipLHC/AdvMuFilterPoint.h
Original file line number Diff line number Diff line change
Expand Up @@ -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)
};

Expand Down
10 changes: 10 additions & 0 deletions shipLHC/AdvTargetHit.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,16 @@ AdvTargetHit::AdvTargetHit(Int_t detID, const std::vector<AdvTargetPoint*>& V)
LOG(DEBUG) << "signal created";
}

AdvTargetHit::AdvTargetHit(Int_t detID, const std::vector<const AdvTargetPoint*>& V)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as for the constructor in AdvMuFilter

: 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
{
Expand Down
17 changes: 9 additions & 8 deletions shipLHC/AdvTargetHit.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ class AdvTargetHit : public SndlhcHit

// Constructor from AdvTargetPoint
AdvTargetHit(Int_t detID, const std::vector<AdvTargetPoint*>&);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same as for the constructor in AdvMuFilter

AdvTargetHit(Int_t detID, const std::vector<const AdvTargetPoint*>&);

/** Destructor **/
~AdvTargetHit() = default;
Expand All @@ -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
Expand Down
19 changes: 8 additions & 11 deletions shipLHC/AdvTargetPoint.h
Original file line number Diff line number Diff line change
Expand Up @@ -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); }

Expand All @@ -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)
};

Expand Down
Loading