diff --git a/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/FT3BaseParam.h b/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/FT3BaseParam.h index b286aa068611c..7160067f075f7 100644 --- a/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/FT3BaseParam.h +++ b/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/FT3BaseParam.h @@ -42,9 +42,6 @@ struct FT3BaseParam : public o2::conf::ConfigurableParamHelper { Float_t etaOut = 1.5; Float_t Layerx2X0 = 0.01; - // FT3Geometry::External file - std::string configFile = ""; // Overrides geoModel parameter when provided - O2ParamDef(FT3BaseParam, "FT3Base"); }; diff --git a/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/GeometryTGeo.h index bb22af7ad7f9d..3c78850dffb55 100644 --- a/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/FT3/base/include/FT3Base/GeometryTGeo.h @@ -101,9 +101,6 @@ class GeometryTGeo : public o2::itsmft::GeometryTGeo static const char* composeSymNameSensor(Int_t d, Int_t lr); protected: - static constexpr int MAXLAYERS = 15; ///< max number of active layers - - Int_t mNumberOfLayers; ///< number of layers static std::string sInnerVolumeName; ///< Mother inner volume name static std::string sVolumeName; ///< Mother volume name static std::string sLayerName; ///< Layer name diff --git a/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/Detector.h b/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/Detector.h index a68f8cf7788b6..8bc4b7f634d7c 100644 --- a/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/Detector.h +++ b/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/Detector.h @@ -25,7 +25,6 @@ #include "TGeoManager.h" // for gGeoManager, TGeoManager (ptr only) #include "TLorentzVector.h" // for TLorentzVector #include "TVector3.h" // for TVector3 -#include "FT3Base/FT3BaseParam.h" class FairVolume; class TGeoVolume; @@ -34,25 +33,10 @@ class TParticle; class TString; -namespace o2 -{ -namespace ft3 +namespace o2::ft3 { class GeometryTGeo; -} -} // namespace o2 -namespace o2 -{ -namespace ft3 -{ -class FT3Layer; -} -} // namespace o2 - -namespace o2 -{ -namespace ft3 -{ +class FT3BaseParam; class FT3Layer; class Detector : public o2::base::DetImpl @@ -108,8 +92,16 @@ class Detector : public o2::base::DetImpl void PostTrack() override { ; } void PreTrack() override { ; } + static constexpr int IdxForwardDisks = 0; + static constexpr int IdxBackwardDisks = 1; /// Returns the number of layers - Int_t getNumberOfLayers() const { return mNumberOfLayers; } + size_t getNumberOfLayers() const + { + if (mLayerName[IdxBackwardDisks].size() != mLayerName[IdxForwardDisks].size()) { + LOG(fatal) << "Number of layers in the two directions are different! Returning 0."; + } + return mLayerName[IdxBackwardDisks].size(); + } void buildBasicFT3(const FT3BaseParam& param); void buildFT3V1(); @@ -117,16 +109,10 @@ class Detector : public o2::base::DetImpl void buildFT3Scoping(); void buildFT3NewVacuumVessel(); void buildFT3ScopingV3(); - void buildFT3FromFile(std::string); - - GeometryTGeo* mGeometryTGeo; //! access to geometry details - - void exportLayout(); protected: std::vector mLayerID; - std::vector> mLayerName; - Int_t mNumberOfLayers; + std::array, 2> mLayerName; // Two sets of layer names, one per direction (forward/backward) private: /// this is transient data about track passing the sensor @@ -154,16 +140,15 @@ class Detector : public o2::base::DetImpl Detector& operator=(const Detector&); - std::vector> mLayers; - bool mIsPipeActivated = true; //! If Alice 3 pipe is present append inner disks to vacuum volume to avoid overlaps + std::array, 2> mLayers; // Two sets of layers, one per direction (forward/backward) + bool mIsPipeActivated = true; //! If Alice 3 pipe is present append inner disks to vacuum volume to avoid overlaps template friend class o2::base::DetImpl; ClassDefOverride(Detector, 1); }; -} // namespace ft3 -} // namespace o2 +} // namespace o2::ft3 #ifdef USESHM namespace o2 diff --git a/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/FT3Layer.h b/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/FT3Layer.h index 44a0ef0f7d8bc..44fd8eb08e444 100644 --- a/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/FT3Layer.h +++ b/Detectors/Upgrades/ALICE3/FT3/simulation/include/FT3Simulation/FT3Layer.h @@ -36,7 +36,7 @@ class FT3Layer : public TObject FT3Layer() = default; // Sample layer constructor - FT3Layer(Int_t layerDirection, Int_t layerNumber, std::string layerName, Float_t z, Float_t rIn, Float_t rOut, Float_t Layerx2X0); + FT3Layer(Int_t layerDirection, Int_t layerNumber, std::string layerName, Float_t z, Float_t rIn, Float_t rOut, Float_t Layerx2X0, bool partOfMiddleLayers); /// Copy constructor FT3Layer(const FT3Layer&) = default; @@ -51,6 +51,7 @@ class FT3Layer : public TObject auto getInnerRadius() const { return mInnerRadius; } auto getOuterRadius() const { return mOuterRadius; } auto getDirection() const { return mDirection; } + bool getIsInMiddleLayer() const { return mIsMiddleLayer; } auto getZ() const { return mZ; } auto getx2X0() const { return mx2X0; } @@ -77,14 +78,15 @@ class FT3Layer : public TObject static TGeoMedium* medFoam; private: - Int_t mLayerNumber = -1; ///< Current layer number - Int_t mDirection; ///< Layer direction 0=Forward 1 = Backward - std::string mLayerName; ///< Current layer name - Double_t mInnerRadius; ///< Inner radius of this layer - Double_t mOuterRadius; ///< Outer radius of this layer - Double_t mZ; ///< Z position of the layer - Double_t mChipThickness; ///< Chip thickness - Double_t mx2X0; ///< Layer material budget x/X0 + Int_t mLayerNumber = -1; ///< Current layer number + Int_t mDirection; ///< Layer direction 0=Forward 1 = Backward + bool mIsMiddleLayer = true; ///< Wether this layer is part of the middle layers + std::string mLayerName; ///< Current layer name + Double_t mInnerRadius; ///< Inner radius of this layer + Double_t mOuterRadius; ///< Outer radius of this layer + Double_t mZ; ///< Z position of the layer + Double_t mChipThickness; ///< Chip thickness + Double_t mx2X0; ///< Layer material budget x/X0 ClassDefOverride(FT3Layer, 0); // ALICE 3 EndCaps geometry }; diff --git a/Detectors/Upgrades/ALICE3/FT3/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/FT3/simulation/src/Detector.cxx index 4b139272834f1..0a93a4061ae44 100644 --- a/Detectors/Upgrades/ALICE3/FT3/simulation/src/Detector.cxx +++ b/Detectors/Upgrades/ALICE3/FT3/simulation/src/Detector.cxx @@ -59,97 +59,6 @@ Detector::Detector() { } -//_________________________________________________________________________________________________ -void Detector::buildFT3FromFile(std::string configFileName) -{ - // Geometry description from file. One line per disk - // z_layer r_in r_out Layerx2X0 - // This simple file reader is not failproof. Do not add empty lines! - - /* - # Sample FT3 configuration - # z_layer r_in r_out Layerx2X0 - -45.3 2.5 9.26 0.0042 - -46.7 2.5 9.26 0.0042 - -48.6 2.5 9.8 0.0042 - -50.0 2.5 9.8 0.0042 - -52.4 2.5 10.43 0.0042 - -53.8 2.5 10.43 0.0042 - -67.7 3.82 13.01 0.0042 - -69.1 3.82 13.01 0.0042 - -76.1 3.92 14.35 0.0042 - -77.5 3.92 14.35 0.0042 - */ - - mLayerName.clear(); - mLayers.clear(); - mLayerID.clear(); - mLayerName.resize(1); - mLayers.resize(1); - - LOG(info) << "Building FT3 Detector: From file"; - LOG(info) << " FT3 detector configuration: " << configFileName; - std::ifstream ifs(configFileName.c_str()); - if (!ifs.good()) { - LOG(fatal) << " Invalid FT3Base.configFile!"; - } - std::string tempstr; - float z_layer, r_in, r_out, Layerx2X0; - char delimiter; - int layerNumber = 0; - while (std::getline(ifs, tempstr)) { - if (tempstr[0] == '#') { - LOG(info) << " Comment: " << tempstr; - continue; - } - std::istringstream iss(tempstr); - iss >> z_layer; - iss >> r_in; - iss >> r_out; - iss >> Layerx2X0; - - int direction = 1; // Forwards - if (z_layer < 0) { - // Backwards - direction = 0; - } - - std::string directionName = std::to_string(direction); - std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); - mLayerName[0].push_back(layerName); - LOG(info) << "Adding Layer " << layerName << " at z = " << z_layer << " ; direction = " << direction << " ; r_in = " << r_in << " ; r_out = " << r_out << " x/X0 = " << Layerx2X0; - auto& thisLayer = mLayers[0].emplace_back(direction, layerNumber, layerName, z_layer, r_in, r_out, Layerx2X0); - layerNumber++; - } - - mNumberOfLayers = layerNumber; - LOG(info) << " Loaded FT3 Detector with " << mNumberOfLayers << " layers"; -} - -//_________________________________________________________________________________________________ -void Detector::exportLayout() -{ - // Export FT3 Layout description to file. - // One line per disk: - // z_layer r_in r_out Layerx2X0 - - std::string configFileName = "FT3_layout.cfg"; - - LOG(info) << "Exporting FT3 Detector layout to " << configFileName; - - std::ofstream fOut(configFileName.c_str(), std::ios::out); - if (!fOut) { - printf("Cannot open file\n"); - return; - } - fOut << "# z_layer r_in r_out Layerx2X0" << std::endl; - for (auto layers_dir : mLayers) { - for (auto layer : layers_dir) { - fOut << layer.getZ() << " " << layer.getInnerRadius() << " " << layer.getOuterRadius() << " " << layer.getx2X0() << std::endl; - } - } -} - //_________________________________________________________________________________________________ void Detector::buildBasicFT3(const FT3BaseParam& param) { @@ -158,28 +67,27 @@ void Detector::buildBasicFT3(const FT3BaseParam& param) LOG(info) << "Building FT3 Detector: Conical Telescope"; - auto z_first = param.z0; - auto z_length = param.zLength; - auto etaIn = param.etaIn; - auto etaOut = param.etaOut; - auto Layerx2X0 = param.Layerx2X0; - mNumberOfLayers = param.nLayers; - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + const int numberOfLayers = param.nLayers; + const auto z_first = param.z0; + const auto z_length = param.zLength; + const auto etaIn = param.etaIn; + const auto etaOut = param.etaOut; + const auto Layerx2X0 = param.Layerx2X0; + mLayerName[IdxBackwardDisks].resize(numberOfLayers); + mLayerName[IdxForwardDisks].resize(numberOfLayers); mLayerID.clear(); - mLayers.resize(2); - for (int direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { - std::string layerName = GeometryTGeo::getFT3LayerPattern() + std::to_string(layerNumber + mNumberOfLayers * direction); + for (int direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { + std::string layerName = GeometryTGeo::getFT3LayerPattern() + std::to_string(layerNumber + numberOfLayers * direction); mLayerName[direction][layerNumber] = layerName; // Adds evenly spaced layers - float layerZ = z_first + (layerNumber * z_length / (mNumberOfLayers - 1)) * std::copysign(1, z_first); - float rIn = std::abs(layerZ * std::tan(2.f * std::atan(std::exp(-etaIn)))); - float rOut = std::abs(layerZ * std::tan(2.f * std::atan(std::exp(-etaOut)))); - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, layerZ, rIn, rOut, Layerx2X0); + const float layerZ = z_first + (layerNumber * z_length / numberOfLayers) * std::copysign(1, z_first); + const float rIn = std::abs(layerZ * std::tan(2.f * std::atan(std::exp(-etaIn)))); + const float rOut = std::abs(layerZ * std::tan(2.f * std::atan(std::exp(-etaOut)))); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, layerZ, rIn, rOut, Layerx2X0, isMiddleLayer); } } } @@ -192,10 +100,10 @@ void Detector::buildFT3V1() LOG(info) << "Building FT3 Detector: V1"; - mNumberOfLayers = 10; - float sensorThickness = 30.e-4; - float layersx2X0 = 1.e-2; - std::vector> layersConfig{ + const int numberOfLayers = 10; + const float sensorThickness = 30.e-4; + const float layersx2X0 = 1.e-2; + const std::vector> layersConfig{ {26., .5, 3., 0.1f * layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} {30., .5, 3., 0.1f * layersx2X0}, {34., .5, 3., 0.1f * layersx2X0}, @@ -207,14 +115,12 @@ void Detector::buildFT3V1() {220., 3.5, 80.f, layersx2X0}, {279., 3.5, 80.f, layersx2X0}}; - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + mLayerName[IdxBackwardDisks].resize(numberOfLayers); + mLayerName[IdxForwardDisks].resize(numberOfLayers); mLayerID.clear(); - mLayers.resize(2); - for (auto direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { + for (auto direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { std::string directionName = std::to_string(direction); std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); mLayerName[direction][layerNumber] = layerName; @@ -226,7 +132,8 @@ void Detector::buildFT3V1() LOG(info) << "Adding Layer " << layerName << " at z = " << z; // Add layers - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0, isMiddleLayer); } } } @@ -239,10 +146,10 @@ void Detector::buildFT3V3b() LOG(info) << "Building FT3 Detector: V3b"; - mNumberOfLayers = 12; + const int numberOfLayers = 12; float sensorThickness = 30.e-4; float layersx2X0 = 1.e-2; - std::vector> layersConfig{ + std::vector> layersConfig{ {26., .5, 3., 0.1f * layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} {30., .5, 3., 0.1f * layersx2X0}, {34., .5, 3., 0.1f * layersx2X0}, @@ -256,14 +163,12 @@ void Detector::buildFT3V3b() {340., 12.5, 80.f, layersx2X0}, {400., 14.7, 80.f, layersx2X0}}; - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + mLayerName[IdxBackwardDisks].resize(numberOfLayers); + mLayerName[IdxForwardDisks].resize(numberOfLayers); mLayerID.clear(); - mLayers.resize(2); - for (auto direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { + for (auto direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { std::string directionName = std::to_string(direction); std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); mLayerName[direction][layerNumber] = layerName; @@ -275,7 +180,8 @@ void Detector::buildFT3V3b() LOG(info) << "Adding Layer " << layerName << " at z = " << z; // Add layers - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0, isMiddleLayer); } } } @@ -291,10 +197,10 @@ void Detector::buildFT3NewVacuumVessel() LOG(info) << "Building FT3 Detector: After Upgrade Days March 2024 version"; - mNumberOfLayers = 9; - float sensorThickness = 30.e-4; - float layersx2X0 = 1.e-2; - std::vector> layersConfigCSide{ + const int numberOfLayers = 9; + const float sensorThickness = 30.e-4; + const float layersx2X0 = 1.e-2; + const std::vector> layersConfigCSide{ {77., 7.0, 35., layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} {100., 7.0, 35., layersx2X0}, {122., 7.0, 35., layersx2X0}, @@ -305,7 +211,7 @@ void Detector::buildFT3NewVacuumVessel() {300., 7.0, 68.f, layersx2X0}, {350., 7.0, 68.f, layersx2X0}}; - std::vector> layersConfigASide{ + const std::vector> layersConfigASide{ {77., 5.0, 35., layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} {100., 5.0, 35., layersx2X0}, {122., 5.0, 35., layersx2X0}, @@ -316,14 +222,12 @@ void Detector::buildFT3NewVacuumVessel() {300., 5.0, 68.f, layersx2X0}, {350., 5.0, 68.f, layersx2X0}}; - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + mLayerName[IdxBackwardDisks].resize(numberOfLayers); + mLayerName[IdxForwardDisks].resize(numberOfLayers); mLayerID.clear(); - mLayers.resize(2); - for (auto direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { + for (auto direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { std::string directionName = std::to_string(direction); std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); mLayerName[direction][layerNumber] = layerName; @@ -342,7 +246,8 @@ void Detector::buildFT3NewVacuumVessel() LOG(info) << "Adding Layer " << layerName << " at z = " << z; // Add layers - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0, isMiddleLayer); } } } @@ -356,52 +261,45 @@ void Detector::buildFT3ScopingV3() LOG(info) << "Building FT3 Detector: v3 scoping version"; - mNumberOfLayers = 6; - float sensorThickness = 30.e-4; - float layersx2X0 = 1.e-2; - std::vector> layersConfigCSide{ - {77., 10.0, 35., layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} - {100., 10.0, 35., layersx2X0}, - {122., 10.0, 35., layersx2X0}, - {150., 20.0, 68.f, layersx2X0}, - {180., 20.0, 68.f, layersx2X0}, - {220., 20.0, 68.f, layersx2X0}}; - - std::vector> layersConfigASide{ - {77., 10.0, 35., layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} - {100., 10.0, 35., layersx2X0}, - {122., 10.0, 35., layersx2X0}, - {150., 20.0, 68.f, layersx2X0}, - {180., 20.0, 68.f, layersx2X0}, - {220., 20.0, 68.f, layersx2X0}}; - - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + const int numberOfLayers = 6; + const float sensorThickness = 30.e-4; + const float layersx2X0 = 1.e-2; + using LayerConfig = std::array; // {z_layer, r_in, r_out, Layerx2X0} + const std::array layersConfigCSide{LayerConfig{77., 10.0, 35., layersx2X0}, + LayerConfig{100., 10.0, 35., layersx2X0}, + LayerConfig{122., 10.0, 35., layersx2X0}, + LayerConfig{150., 20.0, 68.f, layersx2X0}, + LayerConfig{180., 20.0, 68.f, layersx2X0}, + LayerConfig{220., 20.0, 68.f, layersx2X0}}; + + const std::array layersConfigASide{LayerConfig{77., 10.0, 35., layersx2X0}, + LayerConfig{100., 10.0, 35., layersx2X0}, + LayerConfig{122., 10.0, 35., layersx2X0}, + LayerConfig{150., 20.0, 68.f, layersx2X0}, + LayerConfig{180., 20.0, 68.f, layersx2X0}, + LayerConfig{220., 20.0, 68.f, layersx2X0}}; + const std::array enabled{true, true, true, true, true, true}; // To enable or disable layers for debug purpose + mLayerID.clear(); - mLayers.resize(2); - for (auto direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { - std::string directionName = std::to_string(direction); - std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); - mLayerName[direction][layerNumber] = layerName; - float z, rIn, rOut, x0; - if (direction == 0) { // C-Side - z = layersConfigCSide[layerNumber][0]; - rIn = layersConfigCSide[layerNumber][1]; - rOut = layersConfigCSide[layerNumber][2]; - x0 = layersConfigCSide[layerNumber][3]; - } else if (direction == 1) { // A-Side - z = layersConfigASide[layerNumber][0]; - rIn = layersConfigASide[layerNumber][1]; - rOut = layersConfigASide[layerNumber][2]; - x0 = layersConfigASide[layerNumber][3]; + for (int direction : {IdxBackwardDisks, IdxForwardDisks}) { + mLayerName[direction].clear(); + const std::array& layerConfig = (direction == IdxBackwardDisks) ? layersConfigCSide : layersConfigASide; + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { + if (!enabled[layerNumber]) { + continue; } - - LOG(info) << "Adding Layer " << layerName << " at z = " << z; + const std::string directionName = std::to_string(direction); + const std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); + mLayerName[direction].push_back(layerName.c_str()); + const float z = layerConfig[layerNumber][0]; + const float rIn = layerConfig[layerNumber][1]; + const float rOut = layerConfig[layerNumber][2]; + const float x0 = layerConfig[layerNumber][3]; + LOG(info) << "buildFT3ScopingV3 -> Adding Layer " << layerNumber << "/" << numberOfLayers << " " << layerName << " at z = " << z; // Add layers - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0, isMiddleLayer); } } } @@ -413,10 +311,10 @@ void Detector::buildFT3Scoping() LOG(info) << "Building FT3 Detector: Scoping document version"; - mNumberOfLayers = 12; - float sensorThickness = 30.e-4; - float layersx2X0 = 1.e-2; - std::vector> layersConfig{ + const int numberOfLayers = 12; + const float sensorThickness = 30.e-4; + const float layersx2X0 = 1.e-2; + const std::vector> layersConfig{ {26., .5, 2.5, 0.1f * layersx2X0}, // {z_layer, r_in, r_out, Layerx2X0} {30., .5, 2.5, 0.1f * layersx2X0}, {34., .5, 2.5, 0.1f * layersx2X0}, @@ -430,26 +328,24 @@ void Detector::buildFT3Scoping() {300., 5.0, 68.f, layersx2X0}, {350., 5.0, 68.f, layersx2X0}}; - mLayerName.resize(2); - mLayerName[0].resize(mNumberOfLayers); - mLayerName[1].resize(mNumberOfLayers); + mLayerName[IdxBackwardDisks].resize(numberOfLayers); + mLayerName[IdxForwardDisks].resize(numberOfLayers); mLayerID.clear(); - mLayers.resize(2); - for (auto direction : {0, 1}) { - for (int layerNumber = 0; layerNumber < mNumberOfLayers; layerNumber++) { + for (auto direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int layerNumber = 0; layerNumber < numberOfLayers; layerNumber++) { std::string directionName = std::to_string(direction); std::string layerName = GeometryTGeo::getFT3LayerPattern() + directionName + std::string("_") + std::to_string(layerNumber); mLayerName[direction][layerNumber] = layerName; auto& z = layersConfig[layerNumber][0]; - auto& rIn = layersConfig[layerNumber][1]; auto& rOut = layersConfig[layerNumber][2]; auto& x0 = layersConfig[layerNumber][3]; LOG(info) << "Adding Layer " << layerName << " at z = " << z; // Add layers - auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0); + const bool isMiddleLayer = layerNumber < 3; + auto& thisLayer = mLayers[direction].emplace_back(direction, layerNumber, layerName, z, rIn, rOut, x0, isMiddleLayer); } } } @@ -464,24 +360,17 @@ Detector::Detector(bool active) // FT3 Base configuration parameters auto& ft3BaseParam = FT3BaseParam::Instance(); - if (ft3BaseParam.configFile != "") { - LOG(info) << "FT3 Geometry configuration file provided. Overriding FT3Base.geoModel configuration."; - buildFT3FromFile(ft3BaseParam.configFile); - - } else { - switch (ft3BaseParam.geoModel) { - case Default: - buildFT3ScopingV3(); // v3 Dec 25 - break; - case Telescope: - buildBasicFT3(ft3BaseParam); // BasicFT3 = Parametrized telescopic detector (equidistant layers) - break; - default: - LOG(fatal) << "Invalid Geometry.\n"; - break; - } + switch (ft3BaseParam.geoModel) { + case Default: + buildFT3ScopingV3(); // v3 Dec 25 + break; + case Telescope: + buildBasicFT3(ft3BaseParam); // BasicFT3 = Parametrized telescopic detector (equidistant layers) + break; + default: + LOG(fatal) << "Invalid Geometry.\n"; + break; } - exportLayout(); } //_________________________________________________________________________________________________ @@ -494,7 +383,6 @@ Detector::Detector(const Detector& rhs) { mLayerID = rhs.mLayerID; mLayerName = rhs.mLayerName; - mNumberOfLayers = rhs.mNumberOfLayers; } //_________________________________________________________________________________________________ @@ -527,7 +415,6 @@ Detector& Detector::operator=(const Detector& rhs) mLayerID = rhs.mLayerID; mLayerName = rhs.mLayerName; - mNumberOfLayers = rhs.mNumberOfLayers; mLayers = rhs.mLayers; mTrackData = rhs.mTrackData; @@ -543,8 +430,6 @@ void Detector::InitializeO2Detector() // Define the list of sensitive volumes LOG(info) << "Initialize FT3 O2Detector"; - mGeometryTGeo = GeometryTGeo::Instance(); - defineSensitiveVolumes(); } @@ -693,12 +578,10 @@ void Detector::ConstructGeometry() void Detector::createGeometry() { - mGeometryTGeo = GeometryTGeo::Instance(); - TGeoVolume* volFT3 = new TGeoVolumeAssembly(GeometryTGeo::getFT3VolPattern()); TGeoVolume* volIFT3 = new TGeoVolumeAssembly(GeometryTGeo::getFT3InnerVolPattern()); - LOG(info) << "GeometryBuilder::buildGeometry volume name = " << GeometryTGeo::getFT3VolPattern(); + LOG(info) << "FT3: createGeometry volume name = " << GeometryTGeo::getFT3VolPattern(); TGeoVolume* vALIC = gGeoManager->GetVolume("barrel"); if (!vALIC) { @@ -710,69 +593,40 @@ void Detector::createGeometry() LOG(info) << "Running simulation with no beam pipe."; } - LOG(debug) << "FT3 createGeometry: " - << Form("gGeoManager name is %s title is %s", gGeoManager->GetName(), gGeoManager->GetTitle()); - - if (mLayers.size() == 2) { // V1 and telescope - if (!A3IPvac) { - for (int direction : {0, 1}) { // Backward layers at mLayers[0]; Forward layers at mLayers[1] - std::string directionString = direction ? "Forward" : "Backward"; - LOG(info) << "Creating FT3 " << directionString << " layers:"; - for (int iLayer = 0; iLayer < mLayers[direction].size(); iLayer++) { - mLayers[direction][iLayer].createLayer(volFT3); - } - } - vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); - } else { // If beampipe is enabled append inner disks to beampipe filling volume, this should be temporary. - for (int direction : {0, 1}) { - std::string directionString = direction ? "Forward" : "Backward"; - LOG(info) << "Creating FT3 " << directionString << " layers:"; - for (int iLayer = 0; iLayer < mLayers[direction].size(); iLayer++) { - if (iLayer < 3) { - mLayers[direction][iLayer].createLayer(volIFT3); - } else { - mLayers[direction][iLayer].createLayer(volFT3); - } - } - } - A3IPvac->AddNode(volIFT3, 2, new TGeoTranslation(0., 0., 0.)); - vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); - } - - for (auto direction : {0, 1}) { - std::string directionString = direction ? "Forward" : "Backward"; - LOG(info) << "Registering FT3 " << directionString << " LayerIDs:"; + // This will need to adapt to the new scheme + if (!A3IPvac) { + for (int direction : {IdxBackwardDisks, IdxForwardDisks}) { // Backward layers at mLayers[0]; Forward layers at mLayers[1] + const std::string directionString = direction ? "Forward" : "Backward"; + LOG(info) << " Creating FT3 without beampipe " << directionString << " layers:"; for (int iLayer = 0; iLayer < mLayers[direction].size(); iLayer++) { - auto layerID = gMC ? TVirtualMC::GetMC()->VolId(Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), direction, iLayer)) : 0; - mLayerID.push_back(layerID); - LOG(info) << " " << directionString << " layer " << iLayer << " LayerID " << layerID; + mLayers[direction][iLayer].createLayer(volFT3); } } - } - - if (mLayers.size() == 1) { // All layers registered at mLayers[0], used when building from file - LOG(info) << "Creating FT3 layers:"; - if (A3IPvac) { - for (int iLayer = 0; iLayer < mLayers[0].size(); iLayer++) { - if (std::abs(mLayers[0][iLayer].getZ()) < 25) { - mLayers[0][iLayer].createLayer(volIFT3); + vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); + } else { // If beampipe is enabled append inner disks to beampipe filling volume, this should be temporary. + for (int direction : {IdxBackwardDisks, IdxForwardDisks}) { + const std::string directionString = direction ? "Forward" : "Backward"; + LOG(info) << " Creating FT3 " << directionString << " layers:"; + for (int iLayer = 0; iLayer < mLayers[direction].size(); iLayer++) { + LOG(info) << " Creating " << directionString << " layer " << iLayer; + if (mLayers[direction][iLayer].getIsInMiddleLayer()) { // ML disks + mLayers[direction][iLayer].createLayer(volIFT3); } else { - mLayers[0][iLayer].createLayer(volFT3); + mLayers[direction][iLayer].createLayer(volFT3); } } - A3IPvac->AddNode(volIFT3, 2, new TGeoTranslation(0., 0., 0.)); - vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); - } else { - for (int iLayer = 0; iLayer < mLayers[0].size(); iLayer++) { - mLayers[0][iLayer].createLayer(volFT3); - } - vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); } - LOG(info) << "Registering FT3 LayerIDs:"; - for (int iLayer = 0; iLayer < mLayers[0].size(); iLayer++) { - auto layerID = gMC ? TVirtualMC::GetMC()->VolId(Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), 0, iLayer)) : 0; + A3IPvac->AddNode(volIFT3, 2, new TGeoTranslation(0., 0., 0.)); + vALIC->AddNode(volFT3, 2, new TGeoTranslation(0., 30., 0.)); + } + + for (auto direction : {IdxBackwardDisks, IdxForwardDisks}) { + std::string directionString = direction ? "Forward" : "Backward"; + LOG(info) << " Registering FT3 " << directionString << " LayerIDs for " << mLayers[direction].size() << " layers:"; + for (int iLayer = 0; iLayer < mLayers[direction].size(); iLayer++) { + auto layerID = gMC ? TVirtualMC::GetMC()->VolId(Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), direction, iLayer)) : 0; mLayerID.push_back(layerID); - LOG(info) << " mLayerID[" << iLayer << "] = " << layerID; + LOG(info) << " " << directionString << " layer " << iLayer << " LayerID " << layerID; } } } @@ -786,40 +640,34 @@ void Detector::defineSensitiveVolumes() TString volumeName; LOG(info) << "Adding FT3 Sensitive Volumes"; - // The names of the FT3 sensitive volumes have the format: FT3Sensor_(0,1)_(0...sNumberLayers-1) - if (mLayers.size() == 2) { - for (int direction : {0, 1}) { - for (int iLayer = 0; iLayer < mNumberOfLayers; iLayer++) { - volumeName = o2::ft3::GeometryTGeo::getFT3SensorPattern() + std::to_string(iLayer); - if (iLayer < 3) { // ML disks - v = geoManager->GetVolume(Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), direction, iLayer)); - AddSensitiveVolume(v); - } else { // OT disks - for (int sensor_count = 0; sensor_count < MAX_SENSORS; ++sensor_count) { - std::string sensor_name_front = "FT3Sensor_front_" + std::to_string(iLayer) + "_" + std::to_string(direction) + "_" + std::to_string(sensor_count); - std::string sensor_name_back = "FT3Sensor_back_" + std::to_string(iLayer) + "_" + std::to_string(direction) + "_" + std::to_string(sensor_count); - v = geoManager->GetVolume(sensor_name_front.c_str()); - if (v) { - AddSensitiveVolume(v); - } - v = geoManager->GetVolume(sensor_name_back.c_str()); - if (v) { - AddSensitiveVolume(v); - } + for (int direction : {IdxBackwardDisks, IdxForwardDisks}) { + for (int iLayer = 0; iLayer < getNumberOfLayers(); iLayer++) { + LOG(info) << "Adding FT3 Sensitive Volume for direction " << direction << " layer " << iLayer << "/" << getNumberOfLayers(); + volumeName = o2::ft3::GeometryTGeo::getFT3SensorPattern() + std::to_string(iLayer); + if (mLayers[direction][iLayer].getIsInMiddleLayer()) { // ML disks + const std::string sensorName = Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), direction, iLayer); + v = geoManager->GetVolume(sensorName.c_str()); + if (!v) { + geoManager->GetListOfVolumes()->ls(); + LOG(fatal) << "Could not find volume " << sensorName << " for direction " << direction << " layer " << iLayer; + } + AddSensitiveVolume(v); + } else { // OT disks + for (int sensor_count = 0; sensor_count < MAX_SENSORS; ++sensor_count) { + std::string sensor_name_front = "FT3sensor_front_" + std::to_string(iLayer) + "_" + std::to_string(direction) + "_" + std::to_string(sensor_count); + std::string sensor_name_back = "FT3sensor_back_" + std::to_string(iLayer) + "_" + std::to_string(direction) + "_" + std::to_string(sensor_count); + v = geoManager->GetVolume(sensor_name_front.c_str()); + if (v) { + AddSensitiveVolume(v); + } + v = geoManager->GetVolume(sensor_name_back.c_str()); + if (v) { + AddSensitiveVolume(v); } } } } } - - if (mLayers.size() == 1) { - for (int iLayer = 0; iLayer < mLayers[0].size(); iLayer++) { - volumeName = o2::ft3::GeometryTGeo::getFT3SensorPattern() + std::to_string(iLayer); - v = geoManager->GetVolume(Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), mLayers[0][iLayer].getDirection(), iLayer)); - LOG(info) << "Adding FT3 Sensitive Volume => " << v->GetName(); - AddSensitiveVolume(v); - } - } } //_________________________________________________________________________________________________ diff --git a/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Layer.cxx b/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Layer.cxx index 97f42eca6143f..1ad4d1aad1eeb 100644 --- a/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Layer.cxx +++ b/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Layer.cxx @@ -16,7 +16,6 @@ #include "FT3Simulation/FT3Layer.h" #include "FT3Base/GeometryTGeo.h" -#include "FT3Simulation/Detector.h" #include // for LOG @@ -55,20 +54,29 @@ TGeoMedium* FT3Layer::waterMed = nullptr; TGeoMaterial* FT3Layer::foamMat = nullptr; TGeoMedium* FT3Layer::medFoam = nullptr; -FT3Layer::FT3Layer(Int_t layerDirection, Int_t layerNumber, std::string layerName, Float_t z, Float_t rIn, Float_t rOut, Float_t Layerx2X0) +FT3Layer::FT3Layer(Int_t layerDirection, Int_t layerNumber, std::string layerName, Float_t z, Float_t rIn, Float_t rOut, Float_t Layerx2X0, bool partOfMiddleLayers) { // Creates a simple parametrized EndCap layer covering the given // pseudorapidity range at the z layer position mDirection = layerDirection; mLayerNumber = layerNumber; + mIsMiddleLayer = partOfMiddleLayers; mLayerName = layerName; mZ = layerDirection ? std::abs(z) : -std::abs(z); mx2X0 = Layerx2X0; mInnerRadius = rIn; mOuterRadius = rOut; - auto Si_X0 = 9.5; + const double Si_X0 = 9.5; mChipThickness = Layerx2X0 * Si_X0; + // Sanity checks + if (std::isnan(mZ)) { + LOG(fatal) << "FT3 Layer " << mLayerNumber << " has z = NaN, which is not a valid number."; + } + if (mZ < 0.001 && mZ > -0.001) { + LOG(fatal) << "FT3 Layer " << mLayerNumber << " has z = " << mZ << " cm, which is very close to 0."; + } + LOG(info) << "Creating FT3 Layer " << mLayerNumber << " ; direction " << mDirection; LOG(info) << " Using silicon X0 = " << Si_X0 << " to emulate layer radiation length."; LOG(info) << " Layer z = " << mZ << " ; R_in = " << mInnerRadius << " ; R_out = " << mOuterRadius << " ; x2X0 = " << mx2X0 << " ; ChipThickness = " << mChipThickness; @@ -110,8 +118,8 @@ void FT3Layer::createSeparationLayer_waterCooling(TGeoVolume* motherVolume, cons FT3Layer::initialize_mat(); - double carbonFiberThickness = 0.01; - double foamSpacingThickness = 0.5; + const double carbonFiberThickness = 0.01; // cm + const double foamSpacingThickness = 0.5; // cm TGeoTube* carbonFiberLayer = new TGeoTube(mInnerRadius, mOuterRadius, carbonFiberThickness / 2); @@ -122,15 +130,15 @@ void FT3Layer::createSeparationLayer_waterCooling(TGeoVolume* motherVolume, cons carbonFiberLayerVol1->SetLineColor(kGray + 2); carbonFiberLayerVol2->SetLineColor(kGray + 2); - double zSeparation = foamSpacingThickness / 2.0 + carbonFiberThickness / 2.0; + const double zSeparation = foamSpacingThickness / 2.0 + carbonFiberThickness / 2.0; motherVolume->AddNode(carbonFiberLayerVol1, 1, new TGeoTranslation(0, 0, mZ - zSeparation)); motherVolume->AddNode(carbonFiberLayerVol2, 1, new TGeoTranslation(0, 0, mZ + zSeparation)); - double pipeOuterRadius = 0.20; - double kaptonThickness = 0.0025; - double pipeInnerRadius = pipeOuterRadius - kaptonThickness; - double pipeMaxLength = mOuterRadius * 2.0; + const double pipeOuterRadius = 0.20; + const double kaptonThickness = 0.0025; + const double pipeInnerRadius = pipeOuterRadius - kaptonThickness; + const double pipeMaxLength = mOuterRadius * 2.0; int name_it = 0; @@ -199,8 +207,8 @@ void FT3Layer::createSeparationLayer(TGeoVolume* motherVolume, const std::string FT3Layer::initialize_mat(); - double carbonFiberThickness = 0.01; - double foamSpacingThickness = 1.0; + constexpr double carbonFiberThickness = 0.01; // cm + constexpr double foamSpacingThickness = 1.0; // cm TGeoTube* carbonFiberLayer = new TGeoTube(mInnerRadius, mOuterRadius, carbonFiberThickness / 2); TGeoTube* foamLayer = new TGeoTube(mInnerRadius, mOuterRadius, foamSpacingThickness / 2); @@ -215,16 +223,19 @@ void FT3Layer::createSeparationLayer(TGeoVolume* motherVolume, const std::string foamLayerVol->SetFillColorAlpha(kBlack, 1.0); carbonFiberLayerVol2->SetLineColor(kGray + 2); - double zSeparation = foamSpacingThickness / 2.0 + carbonFiberThickness / 2.0; + const double zSeparation = foamSpacingThickness / 2.0 + carbonFiberThickness / 2.0; - motherVolume->AddNode(carbonFiberLayerVol1, 1, new TGeoTranslation(0, 0, mZ - zSeparation)); - motherVolume->AddNode(foamLayerVol, 1, new TGeoTranslation(0, 0, mZ)); - motherVolume->AddNode(carbonFiberLayerVol2, 1, new TGeoTranslation(0, 0, mZ + zSeparation)); + motherVolume->AddNode(carbonFiberLayerVol1, 1, new TGeoTranslation(0, 0, 0 - zSeparation)); + motherVolume->AddNode(foamLayerVol, 1, new TGeoTranslation(0, 0, 0)); + motherVolume->AddNode(carbonFiberLayerVol2, 1, new TGeoTranslation(0, 0, 0 + zSeparation)); } void FT3Layer::createLayer(TGeoVolume* motherVolume) { - if (mLayerNumber >= 0 && mLayerNumber < 3) { + if (mLayerNumber < 0) { + LOG(fatal) << "Invalid layer number " << mLayerNumber << " for FT3 layer."; + } + if (mIsMiddleLayer) { // ML disks std::string chipName = o2::ft3::GeometryTGeo::getFT3ChipPattern() + std::to_string(mLayerNumber), sensName = Form("%s_%d_%d", GeometryTGeo::getFT3SensorPattern(), mDirection, mLayerNumber); @@ -255,7 +266,7 @@ void FT3Layer::createLayer(TGeoVolume* motherVolume) LOG(info) << "Inserting " << layerVol->GetName() << " inside " << motherVolume->GetName(); motherVolume->AddNode(layerVol, 1, FwdDiskCombiTrans); - } else if (mLayerNumber >= 3) { + } else { // OT disks FT3Module module; @@ -264,11 +275,23 @@ void FT3Layer::createLayer(TGeoVolume* motherVolume) std::string backLayerName = o2::ft3::GeometryTGeo::getFT3LayerPattern() + std::to_string(mDirection) + std::to_string(mLayerNumber) + "_Back"; std::string separationLayerName = "FT3SeparationLayer" + std::to_string(mDirection) + std::to_string(mLayerNumber); + TGeoMedium* medAir = gGeoManager->GetMedium("FT3_AIR$"); + TGeoTube* layer = new TGeoTube(mInnerRadius, mOuterRadius, 10 * mChipThickness / 2); + TGeoVolume* layerVol = new TGeoVolume(mLayerName.c_str(), layer, medAir); + layerVol->SetLineColor(kYellow + 2); + // createSeparationLayer_waterCooling(motherVolume, separationLayerName); - createSeparationLayer(motherVolume, separationLayerName); + createSeparationLayer(layerVol, separationLayerName); // create disk faces - module.createModule(mZ, mLayerNumber, mDirection, mInnerRadius, mOuterRadius, 0., "front", "rectangular", motherVolume); - module.createModule(mZ, mLayerNumber, mDirection, mInnerRadius, mOuterRadius, 0., "back", "rectangular", motherVolume); + module.createModule(0, mLayerNumber, mDirection, mInnerRadius, mOuterRadius, 0., "front", "rectangular", layerVol); + module.createModule(0, mLayerNumber, mDirection, mInnerRadius, mOuterRadius, 0., "back", "rectangular", layerVol); + + // Finally put everything in the mother volume + auto* FwdDiskRotation = new TGeoRotation("FwdDiskRotation", 0, 0, 180); + auto* FwdDiskCombiTrans = new TGeoCombiTrans(0, 0, mZ, FwdDiskRotation); + + LOG(info) << "Inserting " << layerVol->GetName() << " inside " << motherVolume->GetName(); + motherVolume->AddNode(layerVol, 1, FwdDiskCombiTrans); } } diff --git a/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Module.cxx b/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Module.cxx index 9318554837706..20a481cb36046 100644 --- a/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Module.cxx +++ b/Detectors/Upgrades/ALICE3/FT3/simulation/src/FT3Module.cxx @@ -199,7 +199,7 @@ void FT3Module::create_layout(double mZ, int layerNumber, int direction, double bottom_y_pos_value = 3.5; bottom_y_neg_value = -3.5; } else { - std::cout << "Different config - to determine offsets needed." << std::endl; + LOG(warning) << "Different config - to determine offsets needed for " << "Rin = " << Rin << " ; sensor_height = " << sensor_height << " ; sensor_width = " << sensor_width << " layer " << layerNumber; x_condition_min = -Rin; x_condition_max = Rin; adjust_bottom_y_pos = false; diff --git a/Detectors/Upgrades/ALICE3/IOTOF/README.md b/Detectors/Upgrades/ALICE3/IOTOF/README.md index 044798076b485..fba4d12252af6 100644 --- a/Detectors/Upgrades/ALICE3/IOTOF/README.md +++ b/Detectors/Upgrades/ALICE3/IOTOF/README.md @@ -14,16 +14,16 @@ Configurables for various sub-detectors are presented in the following Table: [link to definitions](./base/include/IOTOFBase/IOTOFBaseParam.h) -| Options | Choices | Comments | -| ----------------------------- | ---------------------------------------------------------------- | ------------------------------------------- | -| `IOTOFBase.enableInnerTOF` | `true` (default), `false` | Enable inner TOF barrel layer | -| `IOTOFBase.enableOuterTOF` | `true` (default), `false` | Enable outer TOF barrel layer | -| `IOTOFBase.enableForwardTOF` | `true` (default), `false` | Enable forward TOF endcap | -| `IOTOFBase.enableBackwardTOF` | `true` (default), `false` | Enable backward TOF endcap | -| `IOTOFBase.segmentedInnerTOF` | `false` (default), `true` | Use segmented geometry for inner TOF | -| `IOTOFBase.segmentedOuterTOF` | `false` (default), `true` | Use segmented geometry for outer TOF | -| `IOTOFBase.detectorPattern` | ` ` (default), `v3b`, `v3b1a`, `v3b1b`, `v3b2a`, `v3b2b`, `v3b3` | Optional layout pattern | -| ----------------------------- | ------------------------- | ------------------------------------------- | +| Options | Choices | Comments | +| ----------------------------- | ---------------------------------------------------------------- | ---------------------------------------------- | +| `IOTOFBase.enableInnerTOF` | `true` (default), `false` | Enable inner TOF barrel layer | +| `IOTOFBase.enableOuterTOF` | `true` (default), `false` | Enable outer TOF barrel layer | +| `IOTOFBase.enableForwardTOF` | `true` (default), `false` | Enable forward TOF endcap | +| `IOTOFBase.enableBackwardTOF` | `true` (default), `false` | Enable backward TOF endcap | +| `IOTOFBase.segmentedInnerTOF` | `false` (default), `true` | Use segmented geometry for inner TOF | +| `IOTOFBase.segmentedOuterTOF` | `false` (default), `true` | Use segmented geometry for outer TOF | +| `IOTOFBase.detectorPattern` | ` ` (default), `v3b`, `v3b1a`, `v3b1b`, `v3b2a`, `v3b2b`, `v3b3` | Optional layout pattern | +| `IOTOFBase.x2x0` | `0.02` (default) | Chip thickness in fractions of the rad. lenght | For example, a geometry with fully cylindrical tracker barrel (for all layers in VD, ML and OT) can be obtained by diff --git a/Detectors/Upgrades/ALICE3/TRK/README.md b/Detectors/Upgrades/ALICE3/TRK/README.md index 8b3a7984bb233..5b94640cab406 100644 --- a/Detectors/Upgrades/ALICE3/TRK/README.md +++ b/Detectors/Upgrades/ALICE3/TRK/README.md @@ -12,11 +12,11 @@ This is top page for the TRK detector documentation. Configurables for various sub-detectors are presented in the following Table: -| Subsystem | Available options | Comments | -| ------------------ | ------------------------------------------------------- | ---------------------------------------------------------------- | -| `TRKBase.layoutVD` | `kIRIS4` (default), `kIRISFullCyl`, `kIRIS5`, `kIRIS4a` | [link to definitions](./base/include/TRKBase/TRKBaseParam.h) | -| `TRKBase.layoutML` | `kCylinder`, `kTurboStaves` (default), `kStaggered` | | -| `TRKBase.layoutOT` | `kCylinder`, `kTurboStaves`, `kStaggered` (default) | | +| Subsystem | Available options | Comments | +| ------------------ | -------------------------------------------------------------- | ------------------------------------------------------------ | +| `TRKBase.layoutVD` | `kIRIS4` (default), `kIRISFullCyl`, `kIRIS5`, `kIRIS4a` | [link to definitions](./base/include/TRKBase/TRKBaseParam.h) | +| `TRKBase.layoutML` | `kCylinder`, `kTurboStaves` (default), `kStaggered` | | +| `TRKBase.layoutOT` | `kCylinder`, `kTurboStaves`, `kStaggered` (default) | | For example, a geometry with fully cylindrical tracker barrel (for all layers in VD, ML and OT) can be obtained by ```bash diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt index b9866c7d6aa4d..81a75e209124a 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/CMakeLists.txt @@ -9,10 +9,15 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. +if(Acts_FOUND) + set(actsTarget Acts::Core) +endif() + o2_add_library(TRKReconstruction TARGETVARNAME targetName SOURCES src/TimeFrame.cxx src/Clusterer.cxx + $<$:src/ClustererACTS.cxx> PUBLIC_LINK_LIBRARIES O2::ITStracking O2::GPUCommon @@ -27,11 +32,22 @@ o2_add_library(TRKReconstruction O2::DataFormatsITS O2::TRKSimulation nlohmann_json::nlohmann_json + ${actsTarget} PRIVATE_LINK_LIBRARIES O2::Steer TBB::tbb) +if(Acts_FOUND) + target_compile_definitions(${targetName} PUBLIC O2_WITH_ACTS) +endif() + +set(dictHeaders include/TRKReconstruction/TimeFrame.h + include/TRKReconstruction/Clusterer.h) + +if(Acts_FOUND) + list(APPEND dictHeaders include/TRKReconstruction/ClustererACTS.h) +endif() + o2_target_root_dictionary(TRKReconstruction - HEADERS include/TRKReconstruction/TimeFrame.h - include/TRKReconstruction/Clusterer.h + HEADERS ${dictHeaders} LINKDEF src/TRKReconstructionLinkDef.h) diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h index abddafa312fb9..70518b2ace593 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/Clusterer.h @@ -161,17 +161,17 @@ class Clusterer }; //---------------------------------------------- - void process(gsl::span digits, - gsl::span digitROFs, - std::vector& clusters, - std::vector& patterns, - std::vector& clusterROFs, - const ConstDigitTruth* digitLabels = nullptr, - ClusterTruth* clusterLabels = nullptr, - gsl::span digMC2ROFs = {}, - std::vector* clusterMC2ROFs = nullptr); - - private: + virtual void process(gsl::span digits, + gsl::span digitROFs, + std::vector& clusters, + std::vector& patterns, + std::vector& clusterROFs, + const ConstDigitTruth* digitLabels = nullptr, + ClusterTruth* clusterLabels = nullptr, + gsl::span digMC2ROFs = {}, + std::vector* clusterMC2ROFs = nullptr); + + protected: int mNHugeClus = 0; std::unique_ptr mThread; std::vector mSortIdx; ///< reusable per-ROF sort buffer diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h new file mode 100644 index 0000000000000..4111737d17a9f --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/include/TRKReconstruction/ClustererACTS.h @@ -0,0 +1,43 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ClustererACTS.h +/// \brief Definition of the TRK cluster finder + +#ifndef ALICEO2_TRK_CLUSTERERACTS_H +#define ALICEO2_TRK_CLUSTERERACTS_H + +#include "TRKReconstruction/Clusterer.h" + +namespace o2::trk +{ + +class GeometryTGeo; + +class ClustererACTS : public Clusterer +{ + public: + void process(gsl::span digits, + gsl::span digitROFs, + std::vector& clusters, + std::vector& patterns, + std::vector& clusterROFs, + const ConstDigitTruth* digitLabels = nullptr, + ClusterTruth* clusterLabels = nullptr, + gsl::span digMC2ROFs = {}, + std::vector* clusterMC2ROFs = nullptr) override; + + private: +}; + +} // namespace o2::trk + +#endif diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx new file mode 100644 index 0000000000000..711480747ea37 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/ClustererACTS.cxx @@ -0,0 +1,238 @@ +// Copyright 2019-2026 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ClustererACTS.cxx +/// \brief Implementation of the TRK cluster finder with the ACTS + +#include "TRKReconstruction/ClustererACTS.h" +#include "TRKBase/GeometryTGeo.h" +#include + +#include +#include + +using namespace o2::trk; + +struct Cell2D { + Cell2D(int rowv, int colv) : row(rowv), col(colv) {} + int row, col; + Acts::Ccl::Label label{Acts::Ccl::NO_LABEL}; +}; + +int getCellRow(const Cell2D& cell) +{ + return cell.row; +} + +int getCellColumn(const Cell2D& cell) +{ + return cell.col; +} + +bool operator==(const Cell2D& left, const Cell2D& right) +{ + return left.row == right.row && left.col == right.col; +} + +bool cellComp(const Cell2D& left, const Cell2D& right) +{ + return (left.row == right.row) ? left.col < right.col : left.row < right.row; +} + +struct Cluster2D { + std::vector cells; + std::size_t hash{0}; +}; + +void clusterAddCell(Cluster2D& cl, const Cell2D& cell) +{ + cl.cells.push_back(cell); +} + +void hash(Cluster2D& cl) +{ + std::ranges::sort(cl.cells, cellComp); + cl.hash = 0; + // for (const Cell2D& c : cl.cells) { + // boost::hash_combine(cl.hash, c.col); + // } +} + +bool clHashComp(const Cluster2D& left, const Cluster2D& right) +{ + return left.hash < right.hash; +} + +template +void genclusterw(int x, int y, int x0, int y0, int x1, int y1, + std::vector& cells, RNG& rng, double startp = 0.5, + double decayp = 0.9) +{ + std::vector add; + + auto maybe_add = [&](int x_, int y_) { + Cell2D c(x_, y_); + // if (std::uniform_real_distribution()(rng) < startp && + // !rangeContainsValue(cells, c)) { + // cells.push_back(c); + // add.push_back(c); + // } + }; + + // NORTH + if (y < y1) { + maybe_add(x, y + 1); + } + // NORTHEAST + if (x < x1 && y < y1) { + maybe_add(x + 1, y + 1); + } + // EAST + if (x < x1) { + maybe_add(x + 1, y); + } + // SOUTHEAST + if (x < x1 && y > y0) { + maybe_add(x + 1, y - 1); + } + // SOUTH + if (y > y0) { + maybe_add(x, y - 1); + } + // SOUTHWEST + if (x > x0 && y > y0) { + maybe_add(x - 1, y - 1); + } + // WEST + if (x > x0) { + maybe_add(x - 1, y); + } + // NORTHWEST + if (x > x0 && y < y1) { + maybe_add(x - 1, y + 1); + } + + for (Cell2D& c : add) { + genclusterw(c.row, c.col, x0, y0, x1, y1, cells, rng, startp * decayp, + decayp); + } +} + +template +Cluster2D gencluster(int x0, int y0, int x1, int y1, RNG& rng, + double startp = 0.5, double decayp = 0.9) +{ + int x0_ = x0 + 1; + int x1_ = x1 - 1; + int y0_ = y0 + 1; + int y1_ = y1 - 1; + + int x = std::uniform_int_distribution(x0_, x1_)(rng); + int y = std::uniform_int_distribution(y0_, y1_)(rng); + + std::vector cells = {Cell2D(x, y)}; + genclusterw(x, y, x0_, y0_, x1_, y1_, cells, rng, startp, decayp); + + Cluster2D cl; + cl.cells = std::move(cells); + + return cl; +} + +//__________________________________________________ +void ClustererACTS::process(gsl::span digits, + gsl::span digitROFs, + std::vector& clusters, + std::vector& patterns, + std::vector& clusterROFs, + const ConstDigitTruth* digitLabels, + ClusterTruth* clusterLabels, + gsl::span digMC2ROFs, + std::vector* clusterMC2ROFs) +{ + if (!mThread) { + mThread = std::make_unique(this); + } + + auto* geom = o2::trk::GeometryTGeo::Instance(); + + for (size_t iROF = 0; iROF < digitROFs.size(); ++iROF) { + const auto& inROF = digitROFs[iROF]; + const auto outFirst = static_cast(clusters.size()); + const int first = inROF.getFirstEntry(); + const int nEntries = inROF.getNEntries(); + + if (nEntries == 0) { + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), outFirst, 0); + continue; + } + + // template > + // requires(GridDim == 1 || GridDim == 2) + // void createClusters(Acts::Ccl::ClusteringData& data, + // CellCollection& cells, + // ClusterCollection& clusters, + // Connect&& connect = Connect()); + using Cell = Cell2D; + using CellCollection = std::vector; + using Cluster = Cluster2D; + using ClusterCollection = std::vector; + CellCollection cells; + Acts::Ccl::ClusteringData data; + ClusterCollection collection; + + Acts::Ccl::createClusters(data, + cells, + collection, + Acts::Ccl::DefaultConnect(false)); + + // Sort digit indices within this ROF by (chipID, col, row) so we can process + // chip by chip, column by column -- the same ordering the ALPIDE scanner expects. + mSortIdx.resize(nEntries); + std::iota(mSortIdx.begin(), mSortIdx.end(), first); + std::sort(mSortIdx.begin(), mSortIdx.end(), [&digits](int a, int b) { + const auto& da = digits[a]; + const auto& db = digits[b]; + if (da.getChipIndex() != db.getChipIndex()) { + return da.getChipIndex() < db.getChipIndex(); + } + if (da.getColumn() != db.getColumn()) { + return da.getColumn() < db.getColumn(); + } + return da.getRow() < db.getRow(); + }); + + // Process one chip at a time + int sliceStart = 0; + while (sliceStart < nEntries) { + const int chipFirst = sliceStart; + const uint16_t chipID = digits[mSortIdx[sliceStart]].getChipIndex(); + while (sliceStart < nEntries && digits[mSortIdx[sliceStart]].getChipIndex() == chipID) { + ++sliceStart; + } + const int chipN = sliceStart - chipFirst; + + mThread->processChip(digits, chipFirst, chipN, &clusters, &patterns, digitLabels, clusterLabels, geom); + } + + clusterROFs.emplace_back(inROF.getBCData(), inROF.getROFrame(), + outFirst, static_cast(clusters.size()) - outFirst); + } + + if (clusterMC2ROFs && !digMC2ROFs.empty()) { + clusterMC2ROFs->reserve(clusterMC2ROFs->size() + digMC2ROFs.size()); + for (const auto& in : digMC2ROFs) { + clusterMC2ROFs->emplace_back(in.eventRecordID, in.rofRecordID, in.minROF, in.maxROF); + } + } +} diff --git a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h index 4eda22e350852..6ef7321f92e2f 100644 --- a/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h +++ b/Detectors/Upgrades/ALICE3/TRK/reconstruction/src/TRKReconstructionLinkDef.h @@ -17,5 +17,8 @@ #pragma link C++ class o2::trk::TimeFrame < 11> + ; #pragma link C++ class o2::trk::Clusterer + ; +#ifdef O2_WITH_ACTS +#pragma link C++ class o2::trk::ClustererACTS + ; +#endif #endif diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h index bacc1057c7b07..9cfab104ecdf9 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/include/TRKWorkflow/ClustererSpec.h @@ -15,6 +15,9 @@ #include "Framework/DataProcessorSpec.h" #include "Framework/Task.h" #include "TRKReconstruction/Clusterer.h" +#ifdef O2_WITH_ACTS +#include "TRKReconstruction/ClustererACTS.h" +#endif namespace o2::trk { @@ -29,7 +32,13 @@ class ClustererDPL : public o2::framework::Task private: bool mUseMC = true; int mNThreads = 1; +#ifdef O2_WITH_ACTS + bool mUseACTS = false; +#endif o2::trk::Clusterer mClusterer; +#ifdef O2_WITH_ACTS + o2::trk::ClustererACTS mClustererACTS; +#endif }; o2::framework::DataProcessorSpec getClustererSpec(bool useMC); diff --git a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx index 8aec63d69206b..8ce3141fc7e27 100644 --- a/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx +++ b/Detectors/Upgrades/ALICE3/TRK/workflow/src/ClustererSpec.cxx @@ -23,6 +23,9 @@ namespace o2::trk void ClustererDPL::init(o2::framework::InitContext& ic) { mNThreads = std::max(1, ic.options().get("nthreads")); +#ifdef O2_WITH_ACTS + mUseACTS = ic.options().get("useACTS"); +#endif } void ClustererDPL::run(o2::framework::ProcessingContext& pc) @@ -48,15 +51,30 @@ void ClustererDPL::run(o2::framework::ProcessingContext& pc) } o2::base::GeometryManager::loadGeometry("o2sim_geometry.root", false, true); - mClusterer.process(digits, - rofs, - clusters, - patterns, - clusterROFs, - mUseMC ? &labels : nullptr, - clusterLabels.get(), - mc2rofs, - mUseMC ? &clusterMC2ROFs : nullptr); +#ifdef O2_WITH_ACTS + if (mUseACTS) { + mClustererACTS.process(digits, + rofs, + clusters, + patterns, + clusterROFs, + mUseMC ? &labels : nullptr, + clusterLabels.get(), + mc2rofs, + mUseMC ? &clusterMC2ROFs : nullptr); + } else +#endif + { + mClusterer.process(digits, + rofs, + clusters, + patterns, + clusterROFs, + mUseMC ? &labels : nullptr, + clusterLabels.get(), + mc2rofs, + mUseMC ? &clusterMC2ROFs : nullptr); + } pc.outputs().snapshot(o2::framework::Output{"TRK", "COMPCLUSTERS", 0}, clusters); pc.outputs().snapshot(o2::framework::Output{"TRK", "PATTERNS", 0}, patterns); diff --git a/Detectors/Upgrades/ALICE3/macros/ALICE3Field.C b/Detectors/Upgrades/ALICE3/macros/ALICE3Field.C index a721a91ed2dcc..96b7dd4bbe858 100644 --- a/Detectors/Upgrades/ALICE3/macros/ALICE3Field.C +++ b/Detectors/Upgrades/ALICE3/macros/ALICE3Field.C @@ -19,8 +19,8 @@ std::function field() double R2; double B1; double B2; - double beamStart = 500.; //[cm] - double tokGauss = 1. / 0.1; // conversion from Tesla to kGauss + const double beamStart = 500.; //[cm] + const double tokGauss = 1. / 0.1; // conversion from Tesla to kGauss bool isMagAbs = true; @@ -63,4 +63,45 @@ std::function field() b[2] = 0.; } }; +} + +void ALICE3Field() +{ + auto fieldFunc = field(); + // RZ plane visualization + TCanvas* cRZ = new TCanvas("cRZ", "Field in RZ plane", 800, 600); + TH2F* hRZ = new TH2F("hRZ", "Magnetic Field B_z in RZ plane;Z [m];R [m]", 100, -10, 10, 100, -5, 5); + hRZ->SetBit(TH1::kNoStats); // disable stats box + for (int i = 1; i <= hRZ->GetNbinsX(); i++) { + const double Z = hRZ->GetXaxis()->GetBinCenter(i); + for (int j = 1; j <= hRZ->GetNbinsY(); j++) { + const double R = hRZ->GetYaxis()->GetBinCenter(j); + const double pos[3] = {R * 100, 0, Z * 100}; // convert to cm + double b[3] = {0, 0, 0}; + fieldFunc(pos, b); + hRZ->SetBinContent(i, j, b[2]); + } + } + + hRZ->Draw("COLZ"); + cRZ->Update(); + + // XY plane visualization + TCanvas* cXY = new TCanvas("cXY", "Field in XY plane", 800, 600); + TH2F* hXY = new TH2F("hXY", "Magnetic Field B_z in XY plane;X [m];Y [m]", 100, -5, 5, 100, -5, 5); + hXY->SetBit(TH1::kNoStats); // disable stats box + + for (int i = 1; i <= hXY->GetNbinsX(); i++) { + const double X = hXY->GetXaxis()->GetBinCenter(i); + for (int j = 1; j <= hXY->GetNbinsY(); j++) { + const double Y = hXY->GetYaxis()->GetBinCenter(j); + const double pos[3] = {X * 100, Y * 100, 0}; // convert to cm + double b[3] = {0, 0, 0}; + fieldFunc(pos, b); + hXY->SetBinContent(i, j, b[2]); + } + } + + hXY->Draw("COLZ"); + cXY->Update(); } \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/macros/ALICE3FieldShortMagnet.C b/Detectors/Upgrades/ALICE3/macros/ALICE3FieldShortMagnet.C new file mode 100644 index 0000000000000..b12adaa0cc530 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/macros/ALICE3FieldShortMagnet.C @@ -0,0 +1,91 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// Author: J. E. Munoz Mendez jesus.munoz@cern.ch + +std::function field() +{ + return [](const double* x, double* b) { + const double Rc = 185.; //[cm] + const double R1 = 220.; //[cm] + const double R2 = 290.; //[cm] + const double B1 = 2.; //[T] + const double B2 = -Rc * Rc / ((R2 * R2 - R1 * R1) * B1); //[T] + const double beamStart = 370.; //[cm] + const double tokGauss = 1. / 0.1; // conversion from Tesla to kGauss + + const bool isMagAbs = true; + + const double r = sqrt(x[0] * x[0] + x[1] * x[1]); + if ((abs(x[2]) <= beamStart) && (r < Rc)) { // We are inside of the central region + b[0] = 0.; + b[1] = 0.; + b[2] = B1 * tokGauss; + } else if ((abs(x[2]) <= beamStart) && (r >= Rc && r < R1)) { // We are in the transition region + b[0] = 0.; + b[1] = 0.; + b[2] = 0.; + } else if ((abs(x[2]) <= beamStart) && (r >= R1 && r < R2)) { // We are within the magnet + b[0] = 0.; + b[1] = 0.; + if (isMagAbs) { + b[2] = B2 * tokGauss; + } else { + b[2] = 0.; + } + } else { // We are outside of the magnet + b[0] = 0.; + b[1] = 0.; + b[2] = 0.; + } + }; +} + +void ALICE3V3Magnet() +{ + auto fieldFunc = field(); + // RZ plane visualization + TCanvas* cRZ = new TCanvas("cRZ", "Field in RZ plane", 800, 600); + TH2F* hRZ = new TH2F("hRZ", "Magnetic Field B_z in RZ plane;Z [m];R [m]", 100, -10, 10, 100, -5, 5); + hRZ->SetBit(TH1::kNoStats); // disable stats box + for (int i = 1; i <= hRZ->GetNbinsX(); i++) { + const double Z = hRZ->GetXaxis()->GetBinCenter(i); + for (int j = 1; j <= hRZ->GetNbinsY(); j++) { + const double R = hRZ->GetYaxis()->GetBinCenter(j); + const double pos[3] = {R * 100, 0, Z * 100}; // convert to cm + double b[3] = {0, 0, 0}; + fieldFunc(pos, b); + hRZ->SetBinContent(i, j, b[2]); + } + } + + hRZ->Draw("COLZ"); + cRZ->Update(); + + // XY plane visualization + TCanvas* cXY = new TCanvas("cXY", "Field in XY plane", 800, 600); + TH2F* hXY = new TH2F("hXY", "Magnetic Field B_z in XY plane;X [m];Y [m]", 100, -5, 5, 100, -5, 5); + hXY->SetBit(TH1::kNoStats); // disable stats box + + for (int i = 1; i <= hXY->GetNbinsX(); i++) { + const double X = hXY->GetXaxis()->GetBinCenter(i); + for (int j = 1; j <= hXY->GetNbinsY(); j++) { + const double Y = hXY->GetYaxis()->GetBinCenter(j); + const double pos[3] = {X * 100, Y * 100, 0}; // convert to cm + double b[3] = {0, 0, 0}; + fieldFunc(pos, b); + hXY->SetBinContent(i, j, b[2]); + } + } + + hXY->Draw("COLZ"); + cXY->Update(); +} \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/macros/CMakeLists.txt b/Detectors/Upgrades/ALICE3/macros/CMakeLists.txt index b31687cc85c0e..0a4f9031355a2 100644 --- a/Detectors/Upgrades/ALICE3/macros/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/macros/CMakeLists.txt @@ -13,4 +13,7 @@ o2_add_test_root_macro(scanXX0.C LABELS alice3) o2_add_test_root_macro(plotHits.C - LABELS alice3) \ No newline at end of file + LABELS alice3) + +o2_add_test_root_macro(ALICE3FieldShortMagnet.C + LABELS alice3) diff --git a/dependencies/O2Dependencies.cmake b/dependencies/O2Dependencies.cmake index 8addb87a1a16f..b0946e2955369 100644 --- a/dependencies/O2Dependencies.cmake +++ b/dependencies/O2Dependencies.cmake @@ -113,6 +113,9 @@ set_package_properties(fmt PROPERTIES TYPE REQUIRED) find_package(nlohmann_json) set_package_properties(nlohmann_json PROPERTIES TYPE REQUIRED) +find_package(Acts) +set_package_properties(Acts PROPERTIES TYPE OPTIONAL) + find_package(Boost 1.70 COMPONENTS container thread