diff --git a/Detectors/Upgrades/ALICE3/ECal/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ECal/CMakeLists.txt index 83838a01d13f1..cc0a7b0337619 100644 --- a/Detectors/Upgrades/ALICE3/ECal/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/ECal/CMakeLists.txt @@ -10,4 +10,6 @@ # or submit itself to any jurisdiction. add_subdirectory(base) -add_subdirectory(simulation) \ No newline at end of file +add_subdirectory(simulation) +add_subdirectory(reconstruction) +add_subdirectory(DataFormatsECal) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/CMakeLists.txt new file mode 100644 index 0000000000000..3448d6b31029d --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/CMakeLists.txt @@ -0,0 +1,23 @@ +# 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. + +o2_add_library(DataFormatsECal + SOURCES src/Digit.cxx + SOURCES src/MCLabel.cxx + SOURCES src/Cluster.cxx + PUBLIC_LINK_LIBRARIES O2::CommonDataFormat + O2::SimulationDataFormat + AliceO2::InfoLogger) + +o2_target_root_dictionary(DataFormatsECal + HEADERS include/DataFormatsECal/Digit.h + HEADERS include/DataFormatsECal/MCLabel.h + HEADERS include/DataFormatsECal/Cluster.h) diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Cluster.h b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Cluster.h new file mode 100644 index 0000000000000..4a34ef1679f26 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Cluster.h @@ -0,0 +1,85 @@ +// 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. + +/// \file Cluster.h +/// \brief Definition of ECal cluster class +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_CLUSTER_H +#define ALICEO2_ECAL_CLUSTER_H +#include +#include +#include +#include + +namespace o2 +{ +namespace ecal +{ +class Cluster +{ + public: + Cluster() = default; + Cluster(const Cluster& clu) = default; + ~Cluster() = default; + + // setters + void addDigit(int digitIndex, int towerId, double energy); + void setNLM(int nMax) { mNLM = nMax; } + void setE(float energy) { mE = energy; } + void setX(float x) { mX = x; } + void setY(float y) { mY = y; } + void setZ(float z) { mZ = z; } + void setChi2(float chi2) { mChi2 = chi2; } + void setEdgeFlag(bool isEdge) { mEdge = isEdge; } + void addMcTrackID(int mcTrackID, float energy) { mMcTrackEnergy[mcTrackID] += energy; } + + // getters + const std::map& getMcTrackEnergy() { return mMcTrackEnergy; } + int getMultiplicity() const { return mDigitIndex.size(); } + int getDigitIndex(int i) const { return mDigitIndex[i]; } + int getDigitTowerId(int i) const { return mDigitTowerId[i]; } + float getDigitEnergy(int i) const { return mDigitEnergy[i]; } + float getNLM() const { return mNLM; } + float getTime() const { return mTime; } + float getE() const { return mE; } + float getX() const { return mX; } + float getY() const { return mY; } + float getZ() const { return mZ; } + float getR() const { return std::sqrt(mX * mX + mY * mY); } + float getTheta() const { return std::atan2(getR(), mZ); } + float getEta() const { return -std::log(std::tan(getTheta() / 2.)); } + float getPhi() const { return std::atan2(mY, mX); } + float getChi2() const { return mChi2; } + bool isAtTheEdge() const { return mEdge; } + int getMcTrackID() const; + TLorentzVector getMomentum() const; + + private: + std::vector mDigitIndex; // vector of digit indices in digits vector + std::vector mDigitTowerId; // vector of corresponding digit tower Ids + std::vector mDigitEnergy; // vector of corresponding digit energies + std::map mMcTrackEnergy; // MC track indices and corresponding energies + int mNLM = 0; // number of local maxima in the initial cluster + float mTime = 0; // cluster time + float mE = 0; // cluster energy + float mX = 0; // estimated x-coordinate + float mY = 0; // estimated y-ccordinate + float mZ = 0; // estimated z-ccordinate + float mChi2 = 0; // chi2 wrt EM shape + bool mEdge = 0; // set to true if one of cluster digits is at the chamber edge + ClassDefNV(Cluster, 1); +}; +} // namespace ecal +} // namespace o2 + +#endif // ALICEO2_ECAL_CLUSTER_H diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Digit.h b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Digit.h new file mode 100644 index 0000000000000..cc46a64e2cac0 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/Digit.h @@ -0,0 +1,55 @@ +// 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. + +/// \file Digit.h +/// \brief Definition of ECal digit class +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_DIGIT_H +#define ALICEO2_ECAL_DIGIT_H + +#include +#include + +namespace o2 +{ +namespace ecal +{ +class Digit : public o2::dataformats::TimeStamp +{ + public: + Digit() = default; + Digit(int tower, double amplitudeGeV, double time); + ~Digit() = default; + + // setters + void setTower(int tower) { mTower = tower; } + void setAmplitude(double amplitude) { mAmplitudeGeV = amplitude; } + void setEnergy(double energy) { mAmplitudeGeV = energy; } + void setLabel(int label) { mLabel = label; } + + // getters + int getTower() const { return mTower; } + double getAmplitude() const { return mAmplitudeGeV; } + double getEnergy() const { return mAmplitudeGeV; } + int getLabel() const { return mLabel; } + + private: + double mAmplitudeGeV = 0.; ///< Amplitude (GeV) + int32_t mTower = -1; ///< Tower index (absolute cell ID) + int32_t mLabel = -1; ///< Index of the corresponding entry/entries in the MC label array + ClassDefNV(Digit, 1); +}; + +} // namespace ecal +} // namespace o2 +#endif // ALICEO2_ECAL_DIGIT_H diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/MCLabel.h b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/MCLabel.h new file mode 100644 index 0000000000000..762779977ca53 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/include/DataFormatsECal/MCLabel.h @@ -0,0 +1,41 @@ +// 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. + +/// \file MCLabel.h +/// \brief MCLabel class to store MC truth info for ECal +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_MCLABEL_H +#define ALICEO2_ECAL_MCLABEL_H + +#include + +namespace o2 +{ +namespace ecal +{ +class MCLabel : public o2::MCCompLabel +{ + public: + MCLabel() = default; + MCLabel(int trackID, int eventID, int srcID, bool fake, float edep) : o2::MCCompLabel(trackID, eventID, srcID, fake), mEdep(edep) {} + float getEdep() const { return mEdep; } + + private: + float mEdep = 0; // deposited energy + + ClassDefNV(MCLabel, 1); +}; +} // namespace ecal +} // namespace o2 + +#endif // ALICEO2_ECAL_MCLABEL_H diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Cluster.cxx b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Cluster.cxx new file mode 100644 index 0000000000000..77f7d9219ef6b --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Cluster.cxx @@ -0,0 +1,56 @@ +// 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. + +/// \file Cluster.cxx +/// \brief Implementation of ECal cluster class +/// +/// \author Evgeny Kryshen + +#include +#include +#include +#include + +using namespace o2::ecal; + +ClassImp(Cluster); + +//============================================================================== +void Cluster::addDigit(int digitIndex, int towerId, double energy) +{ + mE += energy; + mDigitIndex.push_back(digitIndex); + mDigitTowerId.push_back(towerId); + mDigitEnergy.push_back(energy); +} + +//============================================================================== +int Cluster::getMcTrackID() const +{ + float maxEnergy = 0; + int maxID = 0; + for (const auto& [mcTrackID, energy] : mMcTrackEnergy) { + if (energy > maxEnergy) { + maxEnergy = energy; + maxID = mcTrackID; + } + } + return maxID; +} + +//============================================================================== +TLorentzVector Cluster::getMomentum() const +{ + double r = std::sqrt(mX * mX + mY * mY + mZ * mZ); + if (r == 0) + return TLorentzVector(); + return TLorentzVector(mE * mX / r, mE * mY / r, mE * mZ / r, mE); +} diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/DataFormatsECalLinkDef.h b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/DataFormatsECalLinkDef.h new file mode 100644 index 0000000000000..5b0190aa10d45 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/DataFormatsECalLinkDef.h @@ -0,0 +1,25 @@ +// 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::ecal::Digit + ; +#pragma link C++ class o2::ecal::MCLabel + ; +#pragma link C++ class o2::ecal::Cluster + ; +#pragma link C++ class std::vector < o2::ecal::Digit> + ; +#pragma link C++ class std::vector < o2::ecal::Cluster> + ; +#include "SimulationDataFormat/MCTruthContainer.h" +#pragma link C++ class o2::dataformats::MCTruthContainer < o2::ecal::MCLabel> + ; +#endif diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Digit.cxx b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Digit.cxx new file mode 100644 index 0000000000000..c339c112c6858 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/Digit.cxx @@ -0,0 +1,24 @@ +// 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. + +/// \file Digit.cxx +/// \brief Implementation of ECal digit class +/// +/// \author Evgeny Kryshen + +#include + +using namespace o2::ecal; + +Digit::Digit(int tower, double amplitudeGeV, double time) + : mTower(tower), mAmplitudeGeV(amplitudeGeV), o2::dataformats::TimeStamp(time) +{ +} diff --git a/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/MCLabel.cxx b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/MCLabel.cxx new file mode 100644 index 0000000000000..4dbd2711f1521 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/DataFormatsECal/src/MCLabel.cxx @@ -0,0 +1,19 @@ +// 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. + +/// \file MCLabel.cxx +/// \brief MCLabel class to store MC truth info for ECal +/// +/// \author Evgeny Kryshen + +#include + +ClassImp(o2::ecal::MCLabel); diff --git a/Detectors/Upgrades/ALICE3/ECal/README.md b/Detectors/Upgrades/ALICE3/ECal/README.md index 288040fbd5fd9..ff58683646409 100644 --- a/Detectors/Upgrades/ALICE3/ECal/README.md +++ b/Detectors/Upgrades/ALICE3/ECal/README.md @@ -1,10 +1,10 @@ # ALICE 3 Electromagnetic Calorimenter -This is top page for the ECL detector documentation. +This is top page for the ECAL detector documentation. \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/base/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ECal/base/CMakeLists.txt index 70017cc051e80..b0e1229662653 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/ECal/base/CMakeLists.txt @@ -10,10 +10,14 @@ # or submit itself to any jurisdiction. o2_add_library(ECalBase - SOURCES src/GeometryTGeo.cxx + SOURCES src/Geometry.cxx + src/GeometryTGeo.cxx src/ECalBaseParam.cxx + src/Hit.cxx PUBLIC_LINK_LIBRARIES O2::DetectorsBase) o2_target_root_dictionary(ECalBase - HEADERS include/ECalBase/GeometryTGeo.h - include/ECalBase/ECalBaseParam.h) \ No newline at end of file + HEADERS include/ECalBase/Geometry.h + include/ECalBase/GeometryTGeo.h + include/ECalBase/ECalBaseParam.h + include/ECalBase/Hit.h) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/ECalBaseParam.h b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/ECalBaseParam.h index b8b7c75e2b7d0..aa0de4119914a 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/ECalBaseParam.h +++ b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/ECalBaseParam.h @@ -9,22 +9,45 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file ECalBaseParam.h +/// \brief Geometry parameters configurable via o2-sim --configKeyValues +/// +/// \author Evgeny Kryshen + #ifndef O2_ECAL_BASEPARAM_H #define O2_ECAL_BASEPARAM_H -#include "CommonUtils/ConfigurableParam.h" -#include "CommonUtils/ConfigurableParamHelper.h" +#include +#include namespace o2 { namespace ecal { struct ECalBaseParam : public o2::conf::ConfigurableParamHelper { - float rMin = 125.0; // cm - float rMax = 155.0; // cm - float zLength = 350.0; // cm - - bool enableFwdEndcap = true; + bool enableFwdEndcap = false; + // general ecal barrel settings + double rMin = 125; // cm + double rMax = 155; // cm + double zLength = 350; // cm + int nSuperModules = 4; + // crystal module specification + int nCrystalModulesZ = 31; + int nCrystalModulesPhi = 96; + double crystalAlphaDeg = 0.4; // degrees + double crystalModuleWidth = 1.9; // cm + double crystalModuleLength = 18; // cm + // sampling module specification + int nSamplingModulesZ = 56; + int nSamplingModulesPhi = 67; + double samplingAlphaDeg = 0.4; // degrees + double samplingModuleWidth = 2.7; // cm + double frontPlateThickness = 1.; // cm + double pbLayerThickness = 0.12; // cm + double scLayerThickness = 0.15; // cm + int nSamplingLayers = 80; + // margin in z between crystal modules and sampling modules + double marginCrystalToSampling = 0.1; // cm O2ParamDef(ECalBaseParam, "ECalBase"); }; @@ -32,4 +55,4 @@ struct ECalBaseParam : public o2::conf::ConfigurableParamHelper { } // namespace ecal } // end namespace o2 -#endif \ No newline at end of file +#endif // O2_ECAL_BASEPARAM_H diff --git a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Geometry.h b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Geometry.h new file mode 100644 index 0000000000000..ecfcb5b7cbad6 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Geometry.h @@ -0,0 +1,99 @@ +// 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. + +/// \file Geometry.h +/// \brief Geometry helper class +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_GEOMETRY_H +#define ALICEO2_ECAL_GEOMETRY_H + +#include +#include +#include + +namespace o2 +{ +namespace ecal +{ +class Geometry +{ + public: + static Geometry& instance() + { + static Geometry sGeom; + return sGeom; + } + + int getNcols() const; + int getNrows() const; + std::pair getSectorChamber(int cellId) const; + std::pair getSectorChamber(int iphi, int iz) const; + + void fillFrontFaceCenterCoordinates(); + int getCellID(int moduleId, int sectorId, bool isCrystal); + void detIdToRelIndex(int cellId, int& chamber, int& sector, int& iphi, int& iz) const; + void detIdToGlobalPosition(int detId, double& x, double& y, double& z); + std::pair globalRowColFromIndex(int cellID) const; + bool isCrystal(int cellID); + int areNeighboursVertex(int detId1, int detId2) const; + + double getTanBeta(int i) { return mTanBeta[i]; } + double getFrontFaceZatMinR(int i) { return mFrontFaceZatMinR[i]; } + double getFrontFaceCenterR(int i) { return mFrontFaceCenterR[i]; } + double getFrontFaceCenterZ(int i) { return mFrontFaceCenterZ[i]; } + double getFrontFaceCenterSamplingPhi(int i) { return mFrontFaceCenterSamplingPhi[i]; } + double getFrontFaceCenterCrystalPhi(int i) { return mFrontFaceCenterCrystalPhi[i]; } + double getFrontFaceCenterTheta(int i) { return mFrontFaceCenterTheta[i]; } + double getRMin() { return mRMin; } + double getCrystalModW() { return mCrystalModW; } + double getSamplingModW() { return mSamplingModW; } + double getCrystalAlpha() { return mCrystalAlpha; } + double getSamplingAlpha() { return mSamplingAlpha; } + double getCrystalDeltaPhi() { return 2 * std::atan(mCrystalModW / 2 / mRMin); } + double getSamplingDeltaPhi() { return 2 * std::atan(mSamplingModW / 2 / mRMin); } + double getCrystalPhiMin(); + double getSamplingPhiMin(); + int getNModulesZ() { return mNModulesZ; } + bool isAtTheEdge(int cellId); + + private: + Geometry(); + Geometry(const Geometry&) = delete; + Geometry& operator=(const Geometry&) = delete; + ~Geometry() = default; + double mRMin{0.}; + int mNSuperModules{0}; + int mNCrystalModulesZ{0}; + int mNSamplingModulesZ{0}; + int mNCrystalModulesPhi{0}; + int mNSamplingModulesPhi{0}; + double mCrystalModW{0.}; + double mSamplingModW{0.}; + double mSamplingAlpha{0.}; + double mCrystalAlpha{0.}; + double mMarginCrystalToSampling{0.}; + int mNModulesZ{0}; + std::vector mFrontFaceZatMinR; + std::vector mFrontFaceCenterR; + std::vector mFrontFaceCenterZ; + std::vector mFrontFaceCenterSamplingPhi; + std::vector mFrontFaceCenterCrystalPhi; + std::vector mFrontFaceCenterTheta; + std::vector mTanBeta; + + ClassDefNV(Geometry, 1); +}; +} // namespace ecal +} // namespace o2 + +#endif // ALICEO2_ECAL_GEOMETRY_H diff --git a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/GeometryTGeo.h b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/GeometryTGeo.h index 1cff6dd7d3313..6975a5378a72f 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/GeometryTGeo.h +++ b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/GeometryTGeo.h @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file GeometryTGeo.h +/// \brief Class containing ECal volume naming patterns +/// +/// \author Evgeny Kryshen + #ifndef ALICEO2_ECAL_GEOMETRYTGEO_H #define ALICEO2_ECAL_GEOMETRYTGEO_H @@ -27,21 +32,25 @@ class GeometryTGeo : public o2::detectors::DetMatrixCache static GeometryTGeo* Instance(); static const char* getECalVolPattern() { return sVolumeName.c_str(); } - static const char* getECalSensorPattern() { return sSensorName.c_str(); } + static const char* getECalSectorPattern() { return sSectorName.c_str(); } + static const char* getECalModulePattern() { return sModuleName.c_str(); } static const char* composeSymNameECal() { return Form("%s_%d", o2::detectors::DetID(o2::detectors::DetID::ECL).getName(), 0); } - static const char* composeSymNameSensor(); // A single sensor for the moment + static const char* composeSymNameSector(int s); + static const char* composeSymNameModule(int s, int m); protected: static std::string sVolumeName; - static std::string sSensorName; + static std::string sSectorName; + static std::string sModuleName; private: static std::unique_ptr sInstance; }; } // namespace ecal } // namespace o2 -#endif + +#endif // ALICEO2_ECAL_GEOMETRYTGEO_H diff --git a/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Hit.h b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Hit.h new file mode 100644 index 0000000000000..006b2df5949e6 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/base/include/ECalBase/Hit.h @@ -0,0 +1,85 @@ +// 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. + +/// \file Hit.h +/// \brief MC hit class to store energy loss per cell and per superparent +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_HIT_H +#define ALICEO2_ECAL_HIT_H + +#include +#include + +namespace o2 +{ +namespace ecal +{ +class Hit : public o2::BasicXYZEHit +{ + public: + /// \brief Default constructor + Hit() = default; + + /// \brief Hit constructor + /// + /// Fully defining information of the ECAL point (position, momentum, energy, track, ...) + /// + /// \param trackID Index of the track, defined as parent track entering the ECAL + /// \param cellID ID of the detector cell + /// \param pos Position vector of the point + /// \param mom Momentum vector for the particle at the point + /// \param tof Time of the hit + /// \param eLoss Energy loss + Hit(int trackID, int cellID, const math_utils::Point3D& pos, + const math_utils::Vector3D& mom, float tof, float eLoss) + : o2::BasicXYZEHit(pos.X(), pos.Y(), pos.Z(), tof, eLoss, trackID, 0), + mPvector(mom), + mCellID(cellID) + { + } + + /// \brief Destructor + ~Hit() = default; + + /// \brief Check whether the points are from the same parent and in the same detector volume + /// \return True if points are the same (origin and detector), false otherwise + bool operator==(const Hit& rhs) const; + + /// \brief Sorting points according to parent particle and detector volume + /// \return True if this point is smaller, false otherwise + bool operator<(const Hit& rhs) const; + + /// \brief Get cell ID + /// \return cell ID + int GetCellID() const { return mCellID; } + + private: + math_utils::Vector3D mPvector; ///< Momentum vector + int32_t mCellID; ///< Cell ID (used instead of short detID) + ClassDefNV(Hit, 1); +}; + +} // namespace ecal +} // namespace o2 + +#ifdef USESHM +namespace std +{ +template <> +class allocator : public o2::utils::ShmAllocator +{ +}; +} // namespace std +#endif + +#endif // ALICEO2_ECAL_HIT_H diff --git a/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseLinkDef.h b/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseLinkDef.h index 3bf7ccd32460c..0f0c0637ce2c1 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseLinkDef.h +++ b/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseLinkDef.h @@ -15,8 +15,11 @@ #pragma link off all classes; #pragma link off all functions; +#pragma link C++ class o2::ecal::Geometry + ; #pragma link C++ class o2::ecal::GeometryTGeo + #pragma link C++ class o2::ecal::ECalBaseParam + ; +#pragma link C++ class o2::ecal::Hit + ; #pragma link C++ class o2::conf::ConfigurableParamHelper < o2::ecal::ECalBaseParam> + ; +#pragma link C++ class std::vector < o2::ecal::Hit> + ; #endif \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseParam.cxx b/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseParam.cxx index 54eb2860526b3..6847f42e26346 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseParam.cxx +++ b/Detectors/Upgrades/ALICE3/ECal/base/src/ECalBaseParam.cxx @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include "ECalBase/ECalBaseParam.h" +/// \file ECalBaseParam.cxx +/// \brief Geometry parameters configurable via o2-sim --configKeyValues +/// +/// \author Evgeny Kryshen -O2ParamImpl(o2::ecal::ECalBaseParam); \ No newline at end of file +#include + +O2ParamImpl(o2::ecal::ECalBaseParam); diff --git a/Detectors/Upgrades/ALICE3/ECal/base/src/Geometry.cxx b/Detectors/Upgrades/ALICE3/ECal/base/src/Geometry.cxx new file mode 100644 index 0000000000000..9483b83f19f49 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/base/src/Geometry.cxx @@ -0,0 +1,264 @@ +// 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. + +/// \file Geometry.cxx +/// \brief Geometry helper class +/// +/// \author Evgeny Kryshen + +#include "TMath.h" +#include +#include +#include +#include +#include "CommonConstants/MathConstants.h" +using namespace o2::ecal; +using o2::constants::math::PIHalf; +using o2::constants::math::TwoPI; + +//============================================================================== +Geometry::Geometry() +{ + auto& pars = ECalBaseParam::Instance(); + pars.updateFromFile("o2sim_configuration.ini", "ECalBase"); + pars.printKeyValues(false, true, true, false); + mRMin = pars.rMin; + mNSuperModules = pars.nSuperModules; + mNCrystalModulesZ = pars.nCrystalModulesZ; + mNSamplingModulesZ = pars.nSamplingModulesZ; + mNCrystalModulesPhi = pars.nCrystalModulesPhi; + mNSamplingModulesPhi = pars.nSamplingModulesPhi; + mCrystalModW = pars.crystalModuleWidth; + mSamplingModW = pars.samplingModuleWidth; + mMarginCrystalToSampling = pars.marginCrystalToSampling; + mCrystalAlpha = pars.crystalAlphaDeg * TMath::DegToRad(); + mSamplingAlpha = pars.samplingAlphaDeg * TMath::DegToRad(); + mNModulesZ = 2 * mNSamplingModulesZ + 2 * mNCrystalModulesZ; + fillFrontFaceCenterCoordinates(); +} + +//============================================================================== +int Geometry::getNcols() const +{ + return mNModulesZ; +} + +//============================================================================== +int Geometry::getNrows() const +{ + return mNSuperModules * (mNCrystalModulesPhi > mNSamplingModulesPhi ? mNCrystalModulesPhi : mNSamplingModulesPhi); +} + +//============================================================================== +double Geometry::getCrystalPhiMin() +{ + double superModuleDeltaPhi = TwoPI / mNSuperModules; + double crystalDeltaPhi = getCrystalDeltaPhi(); + return (superModuleDeltaPhi - crystalDeltaPhi * mNCrystalModulesPhi) / 2.; +} + +//============================================================================== +double Geometry::getSamplingPhiMin() +{ + double superModuleDeltaPhi = TwoPI / mNSuperModules; + double samplingDeltaPhi = getSamplingDeltaPhi(); + return (superModuleDeltaPhi - samplingDeltaPhi * mNSamplingModulesPhi) / 2.; +} + +//============================================================================== +void Geometry::fillFrontFaceCenterCoordinates() +{ + if (mFrontFaceCenterR.size() > 0) + return; + mFrontFaceCenterTheta.resize(mNCrystalModulesZ + mNSamplingModulesZ); + mFrontFaceZatMinR.resize(mNCrystalModulesZ + mNSamplingModulesZ); + mFrontFaceCenterR.resize(mNCrystalModulesZ + mNSamplingModulesZ); + mFrontFaceCenterZ.resize(mNCrystalModulesZ + mNSamplingModulesZ); + mTanBeta.resize(mNCrystalModulesZ + mNSamplingModulesZ); + mFrontFaceCenterSamplingPhi.resize(mNSuperModules * mNSamplingModulesPhi); + mFrontFaceCenterCrystalPhi.resize(mNSuperModules * mNCrystalModulesPhi); + + double superModuleDeltaPhi = TwoPI / mNSuperModules; + double crystalDeltaPhi = getCrystalDeltaPhi(); + double samplingDeltaPhi = getSamplingDeltaPhi(); + double crystalPhiMin = getCrystalPhiMin(); + double samplingPhiMin = getSamplingPhiMin(); + for (int ism = 0; ism < mNSuperModules; ism++) { + // crystal + for (int i = 0; i < mNCrystalModulesPhi; i++) { + double phi0 = superModuleDeltaPhi * ism + crystalPhiMin + crystalDeltaPhi * i; + mFrontFaceCenterCrystalPhi[ism * mNCrystalModulesPhi + i] = phi0; + } + // sampling + for (int i = 0; i < mNSamplingModulesPhi; i++) { + double phi0 = superModuleDeltaPhi * ism + samplingPhiMin + samplingDeltaPhi * i; + mFrontFaceCenterSamplingPhi[ism * mNSamplingModulesPhi + i] = phi0; + } + } + + double theta0 = PIHalf - mCrystalAlpha; + double zAtMinR = mCrystalModW * std::cos(mCrystalAlpha); + + for (int m = 0; m < mNCrystalModulesZ; m++) { + mTanBeta[m] = std::sin(theta0 - mCrystalAlpha) * mCrystalModW / 2 / mRMin; + ROOT::Math::Polar2DVector vMid21(mCrystalModW / 2., PIHalf + theta0); + ROOT::Math::XYPoint pAtMinR(zAtMinR, mRMin); + ROOT::Math::XYPoint pc = pAtMinR + vMid21; + mFrontFaceZatMinR[m] = zAtMinR; + mFrontFaceCenterZ[m] = pc.x(); + mFrontFaceCenterR[m] = pc.y(); + mFrontFaceCenterTheta[m] = theta0; + theta0 -= 2 * mCrystalAlpha; + zAtMinR += mCrystalModW * std::cos(mCrystalAlpha) / std::sin(theta0 + mCrystalAlpha); + } + + theta0 = mFrontFaceCenterTheta[mNCrystalModulesZ - 1] - mCrystalAlpha - mSamplingAlpha; + zAtMinR = mFrontFaceZatMinR[mNCrystalModulesZ - 1]; + zAtMinR += mSamplingModW * std::cos(mSamplingAlpha) / std::sin(theta0 + mSamplingAlpha); + zAtMinR += mMarginCrystalToSampling; + + for (int m = 0; m < mNSamplingModulesZ; m++) { + int i = m + mNCrystalModulesZ; + mTanBeta[i] = std::sin(theta0 - mSamplingAlpha) * mSamplingModW / 2 / mRMin; + ROOT::Math::Polar2DVector vMid21(mSamplingModW / 2., PIHalf + theta0); + ROOT::Math::XYPoint pAtMinR(zAtMinR, mRMin); + ROOT::Math::XYPoint pc = pAtMinR + vMid21; + mFrontFaceZatMinR[i] = zAtMinR; + mFrontFaceCenterZ[i] = pc.x(); + mFrontFaceCenterR[i] = pc.y(); + mFrontFaceCenterTheta[i] = theta0; + theta0 -= 2 * mSamplingAlpha; + zAtMinR += mSamplingModW * std::cos(mSamplingAlpha) / std::sin(theta0 + mSamplingAlpha); + } +} + +int Geometry::getCellID(int moduleId, int sectorId, bool isCrystal) +{ + int cellID = 0; + if (isCrystal) { + if (moduleId % 2 == 0) { // crystal at positive eta + cellID = sectorId * mNModulesZ + moduleId / 2 + mNSamplingModulesZ + mNCrystalModulesZ; + } else { // crystal at negative eta + cellID = sectorId * mNModulesZ - moduleId / 2 + mNSamplingModulesZ + mNCrystalModulesZ - 1; + } + } else { + if (sectorId % 2 == 0) { // sampling at positive eta + cellID = sectorId / 2 * mNModulesZ + moduleId + mNSamplingModulesZ + mNCrystalModulesZ * 2; + } else { // sampling at negative eta + cellID = sectorId / 2 * mNModulesZ - moduleId + mNSamplingModulesZ; + } + } + return cellID; +} + +//============================================================================== +std::pair Geometry::globalRowColFromIndex(int cellID) const +{ + int ip = cellID / mNModulesZ; // row + int iz = cellID % mNModulesZ; // col + return {ip, iz}; +} + +//============================================================================== +bool Geometry::isCrystal(int cellID) +{ + auto [row, col] = globalRowColFromIndex(cellID); + return (col >= mNSamplingModulesZ && col < mNSamplingModulesZ + 2 * mNCrystalModulesZ); +} + +//============================================================================== +std::pair Geometry::getSectorChamber(int cellId) const +{ + int iphi = cellId / mNModulesZ; + int iz = cellId % mNModulesZ; + return getSectorChamber(iphi, iz); +} + +//============================================================================== +std::pair Geometry::getSectorChamber(int iphi, int iz) const +{ + int chamber = iz < mNSamplingModulesZ ? 0 : (iz < mNSamplingModulesZ + 2 * mNCrystalModulesZ ? 1 : 2); + int sector = iphi / (chamber == 1 ? mNCrystalModulesPhi : mNSamplingModulesPhi); + return {sector, chamber}; +} + +//============================================================================== +void Geometry::detIdToRelIndex(int cellId, int& chamber, int& sector, int& iphi, int& iz) const +{ + // 3 chambers - sampling/crystal/sampling + iphi = cellId / mNModulesZ; + iz = cellId % mNModulesZ; + auto pair = getSectorChamber(iphi, iz); + sector = pair.first; + chamber = pair.second; +} + +//============================================================================== +void Geometry::detIdToGlobalPosition(int detId, double& x, double& y, double& z) +{ + int chamber, sector, iphi, iz; + detIdToRelIndex(detId, chamber, sector, iphi, iz); + if (iz < mNSamplingModulesZ + mNCrystalModulesZ) { + z = -mFrontFaceCenterZ[mNSamplingModulesZ + mNCrystalModulesZ - iz - 1]; + } else { + z = +mFrontFaceCenterZ[iz % (mNSamplingModulesZ + mNCrystalModulesZ)]; + } + double phi = chamber == 1 ? mFrontFaceCenterCrystalPhi[iphi] : mFrontFaceCenterSamplingPhi[iphi]; + double r = mFrontFaceCenterR[iz % (mNSamplingModulesZ + mNCrystalModulesZ)]; + x = r * std::cos(phi); + y = r * std::sin(phi); +} + +//============================================================================== +int Geometry::areNeighboursVertex(int detId1, int detId2) const +{ + int ch1, sector1, iphi1, iz1; + int ch2, sector2, iphi2, iz2; + detIdToRelIndex(detId1, ch1, sector1, iphi1, iz1); + detIdToRelIndex(detId2, ch2, sector2, iphi2, iz2); + if (sector1 != sector2 || ch1 != ch2) + return 0; + if (std::abs(iphi1 - iphi2) <= 1 && std::abs(iz1 - iz2) <= 1) + return 1; + return 0; +} + +//============================================================================== +bool Geometry::isAtTheEdge(int cellId) +{ + auto [row, col] = globalRowColFromIndex(cellId); + if (col == 0) + return 1; + if (col == mNSamplingModulesZ) + return 1; + if (col == mNSamplingModulesZ - 1) + return 1; + if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ) + return 1; + if (col == mNSamplingModulesZ + 2 * mNCrystalModulesZ - 1) + return 1; + if (col == mNModulesZ - 1) + return 1; + for (int m = 0; m <= mNSuperModules; m++) { + if (isCrystal(cellId)) { + if (row == m * mNCrystalModulesPhi) + return 1; + if (row == m * mNCrystalModulesPhi - 1) + return 1; + } else { + if (row == m * mNSamplingModulesPhi) + return 1; + if (row == m * mNSamplingModulesPhi - 1) + return 1; + } + } + return 0; +} diff --git a/Detectors/Upgrades/ALICE3/ECal/base/src/GeometryTGeo.cxx b/Detectors/Upgrades/ALICE3/ECal/base/src/GeometryTGeo.cxx index 49f57d8a8c5cc..aca4f5548dc51 100644 --- a/Detectors/Upgrades/ALICE3/ECal/base/src/GeometryTGeo.cxx +++ b/Detectors/Upgrades/ALICE3/ECal/base/src/GeometryTGeo.cxx @@ -9,6 +9,11 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +/// \file GeometryTGeo.cxx +/// \brief Class containing ECal volume naming patterns +/// +/// \author Evgeny Kryshen + #include #include @@ -19,7 +24,8 @@ namespace ecal std::unique_ptr GeometryTGeo::sInstance; std::string GeometryTGeo::sVolumeName = "ECALV"; -std::string GeometryTGeo::sSensorName = "ECALSensor"; +std::string GeometryTGeo::sSectorName = "ECALSector"; +std::string GeometryTGeo::sModuleName = "ECALModule"; GeometryTGeo::GeometryTGeo(bool build, int loadTrans) : DetMatrixCache() { @@ -57,10 +63,15 @@ GeometryTGeo* GeometryTGeo::Instance() return sInstance.get(); } -const char* GeometryTGeo::composeSymNameSensor() +const char* GeometryTGeo::composeSymNameSector(int s) +{ + return Form("%s/%s_%d", composeSymNameECal(), getECalSectorPattern(), s); +} + +const char* GeometryTGeo::composeSymNameModule(int s, int m) { - return Form("%s/%d", composeSymNameECal(), 0); + return Form("%s/%s_%d", composeSymNameSector(s), getECalModulePattern(), m); } } // namespace ecal -} // namespace o2 \ No newline at end of file +} // namespace o2 diff --git a/Detectors/Upgrades/ALICE3/ECal/base/src/Hit.cxx b/Detectors/Upgrades/ALICE3/ECal/base/src/Hit.cxx new file mode 100644 index 0000000000000..ee2034314d2d8 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/base/src/Hit.cxx @@ -0,0 +1,34 @@ +// 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. + +/// \file Hit.cxx +/// \brief MC hit class to store energy loss per cell and per superparent +/// +/// \author Evgeny Kryshen + +#include + +ClassImp(o2::ecal::Hit); + +using namespace o2::ecal; + +bool Hit::operator<(const Hit& rhs) const +{ + if (GetTrackID() != rhs.GetTrackID()) { + return GetTrackID() < rhs.GetTrackID(); + } + return GetCellID() < rhs.GetCellID(); +} + +bool Hit::operator==(const Hit& rhs) const +{ + return (GetCellID() == rhs.GetCellID()) && (GetTrackID() == rhs.GetTrackID()); +} diff --git a/Detectors/Upgrades/ALICE3/ECal/reconstruction/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ECal/reconstruction/CMakeLists.txt new file mode 100644 index 0000000000000..f51a9c067d6b3 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/reconstruction/CMakeLists.txt @@ -0,0 +1,19 @@ +# 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. + +o2_add_library(ECalReconstruction + SOURCES src/Clusterizer.cxx + PUBLIC_LINK_LIBRARIES O2::ECalBase + O2::DataFormatsECal + AliceO2::InfoLogger) + +o2_target_root_dictionary(ECalReconstruction + HEADERS include/ECalReconstruction/Clusterizer.h) diff --git a/Detectors/Upgrades/ALICE3/ECal/reconstruction/include/ECalReconstruction/Clusterizer.h b/Detectors/Upgrades/ALICE3/ECal/reconstruction/include/ECalReconstruction/Clusterizer.h new file mode 100644 index 0000000000000..3bb7cab4b11e3 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/reconstruction/include/ECalReconstruction/Clusterizer.h @@ -0,0 +1,75 @@ +// 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. + +/// \file Clusterizer.h +/// \brief Class for cluster finding and unfolding +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_CLUSTERIZER_H +#define ALICEO2_ECAL_CLUSTERIZER_H + +#include +#include +#include +#include + +using o2::ecal::Cluster; +using o2::ecal::Digit; +class TF1; + +namespace o2 +{ +namespace ecal +{ +class Clusterizer +{ + public: + Clusterizer(bool applyCorrectionZ = 1, bool applyCorrectionE = 1); + ~Clusterizer() = default; + void initialize() {}; + void addDigitToCluster(Cluster& cluster, int row, int col, const gsl::span& digits); + void findClusters(const gsl::span& digits, std::vector& foundClusters, std::vector& unfoldedClusters); + void makeClusters(const gsl::span& digits, std::vector& clusters); + void makeUnfoldings(std::vector& foundClusters, std::vector& unfoldedClusters); + void unfoldOneCluster(Cluster* iniClu, int nMax, int* digitId, float* maxAtEnergy, std::vector& unfoldedClusters); + void evalClusters(std::vector& clusters); + int getNumberOfLocalMax(Cluster& clu, int* maxAt, float* maxAtEnergy); + double showerShape(double dx, double dz, bool isCrystal); + void setLogWeight(double logWeight) { mLogWeight = logWeight; } + void setClusteringThreshold(double threshold) { mClusteringThreshold = threshold; } + void setCrystalDigitThreshold(double threshold) { mCrystalDigitThreshold = threshold; } + void setSamplingDigitThreshold(double threshold) { mSamplingDigitThreshold = threshold; } + + private: + std::vector> mDigitIndices; // 2D map of digit indices used for recursive cluster finding + bool mUnfoldClusters = true; // to perform cluster unfolding + double mCrystalDigitThreshold = 0.040; // minimal energy of crystal digit + double mSamplingDigitThreshold = 0.100; // minimal energy of sampling digit + double mClusteringThreshold = 0.050; // minimal energy of digit to start clustering (GeV) + double mClusteringTimeGate = 1e9; // maximal time difference between digits to be accepted to clusters (in ns) + int mNLMMax = 30; // maximal number of local maxima in unfolding + double mLogWeight = 4.; // cutoff used in log. weight calculation + double mUnfogingEAccuracy = 1.e-4; // accuracy of energy calculation in unfoding prosedure (GeV) + double mUnfogingXZAccuracy = 1.e-2; // accuracy of position calculation in unfolding procedure (cm) + int mNMaxIterations = 100; // maximal number of iterations in unfolding procedure + double mLocalMaximumCut = 0.015; // minimal height of local maximum over neighbours + bool mApplyCorrectionZ = 1; // z-correction + bool mApplyCorrectionE = 1; // energy-correction + TF1* fCrystalShowerShape; //! Crystal shower shape + TF1* fSamplingShowerShape; //! Sampling shower shape + TF1* fCrystalRMS; //! Crystal RMS +}; + +} // namespace ecal +} // namespace o2 + +#endif // ALICEO2_ECAL_CLUSTERIZER_H diff --git a/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/Clusterizer.cxx b/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/Clusterizer.cxx new file mode 100644 index 0000000000000..c84f62b60ec38 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/Clusterizer.cxx @@ -0,0 +1,455 @@ +// 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. + +/// \file Clusterizer.cxx +/// \brief Class for cluster finding and unfolding +/// +/// \author Evgeny Kryshen + +#include +#include +#include +#include +#include +#include +#include + +using namespace o2::ecal; + +//============================================================================== +Clusterizer::Clusterizer(bool applyCorrectionZ, bool applyCorrectionE) +{ + auto& geo = Geometry::instance(); + mDigitIndices.resize(geo.getNrows(), std::vector(geo.getNcols(), -1)); + mApplyCorrectionZ = applyCorrectionZ; + mApplyCorrectionE = applyCorrectionE; + fCrystalShowerShape = new TF1("fCrystal", "x<[1] ? [0]*exp([3]*x+[4]*x*x+[5]*x*x*x) : (x<[2] ? [0]*[6]*exp([7]*x+[8]*x*x) : [0]*[9]*exp([10]*x+[11]*x*x))", 0, 15); + double pc[12]; + pc[0] = 1. / 13.15; + pc[1] = 2.2; + pc[2] = 5; + pc[3] = 4.38969; + pc[4] = -5.15975; + pc[5] = 1.18978; + pc[6] = 1.48726; + pc[7] = -1.54621; + pc[8] = 0.0814617; + pc[9] = 0.0369055; + pc[10] = -0.174372; + pc[11] = -0.0455978; + + fCrystalShowerShape->SetParameters(pc); + + fSamplingShowerShape = new TF1("fSampling", "x<[1] ? [0]*exp([3]*x+[4]*x*x+[5]*x*x*x) : (x<[2] ? [0]*[6]*exp([7]*x+[8]*x*x) : [0]*[9]*exp([10]*x+[11]*x*x))", 0, 15); + double ps[12]; + ps[0] = 1 / 35.6; + ps[1] = 3.2; + ps[2] = 6; + ps[3] = 3.06543; + ps[4] = -2.23235; + ps[5] = 0.325344; + ps[6] = 6.0733; + ps[7] = -1.62713; + ps[8] = 0.0965569; + ps[9] = 0.0765706; + ps[10] = -0.217398; + ps[11] = -0.0204646; + fSamplingShowerShape->SetParameters(ps); + + fCrystalRMS = new TF1("fCrystalRMS", "[0]*x*exp([1]*x+[2]*x*x+[3]*x*x*x)", 0, 2.2); + double p[4]; + p[0] = 1.39814; + p[1] = -6.05426; + p[2] = 6.26678; + p[3] = -1.97092; + fCrystalRMS->SetParameters(p); +} + +//============================================================================== +void Clusterizer::findClusters(const gsl::span& digits, std::vector& foundClusters, std::vector& unfoldedClusters) +{ + foundClusters.clear(); + unfoldedClusters.clear(); + + // Collect list of clusters + makeClusters(digits, foundClusters); + + // Split clusters with several local maxima if necessary + makeUnfoldings(foundClusters, unfoldedClusters); + + // Evaluate cluster position, dispersion etc. + evalClusters(foundClusters); + evalClusters(unfoldedClusters); +} + +//============================================================================== +void Clusterizer::addDigitToCluster(Cluster& cluster, int row, int col, const gsl::span& digits) +{ + auto& geo = Geometry::instance(); + if (row < 0 || row >= geo.getNrows() || col < 0 || col >= geo.getNcols()) + return; + int digitIndex = mDigitIndices[row][col]; + LOGP(debug, " checking row={} and col={} digitIndex={}", row, col, digitIndex); + if (digitIndex < 0) + return; + + const Digit& digit = digits[digitIndex]; + if (cluster.getMultiplicity() > 0) { + // check if new digit is in the same chamber and sector + const Digit& digit2 = digits[cluster.getDigitIndex(0)]; + auto [sector1, ch1] = geo.getSectorChamber(digit.getTower()); + auto [sector2, ch2] = geo.getSectorChamber(digit2.getTower()); + LOGP(debug, " checking if sector and chamber are the same: ({},{}) ({},{})", sector1, ch1, sector2, ch2); + if (sector1 != sector2 || ch1 != ch2) + return; + } + + mDigitIndices[row][col] = -1; + cluster.addDigit(digitIndex, digit.getTower(), digit.getEnergy()); + LOGP(debug, " adding new digit at row={} and col={}", row, col); + addDigitToCluster(cluster, row - 1, col, digits); + addDigitToCluster(cluster, row + 1, col, digits); + addDigitToCluster(cluster, row, col - 1, digits); + addDigitToCluster(cluster, row, col + 1, digits); +} + +//============================================================================== +void Clusterizer::makeClusters(const gsl::span& digits, std::vector& clusters) +{ + // Combine digits into cluster + + int nDigits = digits.size(); + + // reset mDigitIndices + for (auto& rows : mDigitIndices) { + rows.assign(rows.size(), -1); + } + + // fill mDigitIndices + auto& geo = Geometry::instance(); + for (int i = 0; i < nDigits; i++) { + const Digit& digit = digits[i]; + auto [row, col] = geo.globalRowColFromIndex(digit.getTower()); + bool isCrystal = geo.isCrystal(digit.getTower()); + if (isCrystal) { + if (digit.getEnergy() < mCrystalDigitThreshold) + continue; + } else { + if (digit.getEnergy() < mSamplingDigitThreshold) + continue; + } + mDigitIndices[row][col] = i; + } + + // add digit seeds to clusters and recursively add neighbours + for (int i = 0; i < nDigits; i++) { + const Digit& digitSeed = digits[i]; + auto [row, col] = geo.globalRowColFromIndex(digitSeed.getTower()); + if (mDigitIndices[row][col] < 0) + continue; // digit was already added in one of the clusters + if (digitSeed.getEnergy() < mClusteringThreshold) + continue; + LOGP(debug, " starting new cluster at row={} and col={}", row, col); + auto& cluster = clusters.emplace_back(); + addDigitToCluster(cluster, row, col, digits); + } + + LOGP(debug, "made {} clusters from {} digits", clusters.size(), nDigits); +} + +//============================================================================== +void Clusterizer::makeUnfoldings(std::vector& foundClusters, std::vector& unfoldedClusters) +{ + // Split cluster if several local maxima are found + if (!mUnfoldClusters) { + return; + } + + int* maxAt = new int[mNLMMax]; + float* maxAtEnergy = new float[mNLMMax]; + + for (auto& clu : foundClusters) { + int nMax = getNumberOfLocalMax(clu, maxAt, maxAtEnergy); + if (nMax > 1) { + unfoldOneCluster(&clu, nMax, maxAt, maxAtEnergy, unfoldedClusters); + } else { + clu.setNLM(1); + unfoldedClusters.emplace_back(clu); + } + } + delete[] maxAt; + delete[] maxAtEnergy; +} + +//============================================================================== +void Clusterizer::unfoldOneCluster(Cluster* iniClu, int nMax, int* digitId, float* maxAtEnergy, std::vector& unfoldedClusters) +{ + // Based on MpdEmcClusterizerKI::UnfoldOneCluster by D. Peresunko + // Performs the unfolding of a cluster with nMax overlapping showers + // Parameters: iniClu cluster to be unfolded + // nMax number of local maxima found (this is the number of new clusters) + // digitId: index of digits, corresponding to local maxima + // maxAtEnergy: energies of digits, corresponding to local maxima + + // Take initial cluster and calculate local coordinates of digits + // To avoid multiple re-calculation of same parameters + int mult = iniClu->getMultiplicity(); + std::vector x(mult); + std::vector y(mult); + std::vector z(mult); + std::vector e(mult); + std::vector> eInClusters(mult, std::vector(nMax)); + + auto& geo = Geometry::instance(); + bool isCrystal = geo.isCrystal(iniClu->getDigitTowerId(0)); + + for (int idig = 0; idig < mult; idig++) { + e[idig] = iniClu->getDigitEnergy(idig); + geo.detIdToGlobalPosition(iniClu->getDigitTowerId(idig), x[idig], y[idig], z[idig]); + } + + // Coordinates of centers of clusters + std::vector xMax(nMax); + std::vector yMax(nMax); + std::vector zMax(nMax); + std::vector eMax(nMax); + + for (int iclu = 0; iclu < nMax; iclu++) { + xMax[iclu] = x[digitId[iclu]]; + yMax[iclu] = y[digitId[iclu]]; + zMax[iclu] = z[digitId[iclu]]; + eMax[iclu] = e[digitId[iclu]]; + } + + std::vector prop(nMax); // proportion of clusters in the current digit + + // Try to decompose cluster to contributions + int nIterations = 0; + bool insuficientAccuracy = true; + + while (insuficientAccuracy && nIterations < mNMaxIterations) { + // Loop over all digits of parent cluster and split their energies between daughter clusters + // according to shower shape + for (int idig = 0; idig < mult; idig++) { + double eEstimated = 0; + for (int iclu = 0; iclu < nMax; iclu++) { + prop[iclu] = eMax[iclu] * showerShape(std::sqrt((x[idig] - xMax[iclu]) * (x[idig] - xMax[iclu]) + + (y[idig] - yMax[iclu]) * (y[idig] - yMax[iclu])), + z[idig] - zMax[iclu], isCrystal); + eEstimated += prop[iclu]; + } + if (eEstimated == 0.) { // numerical accuracy + continue; + } + // Split energy of digit according to contributions + for (int iclu = 0; iclu < nMax; iclu++) { + eInClusters[idig][iclu] = e[idig] * prop[iclu] / eEstimated; + } + } + // Recalculate parameters of clusters and check relative variation of energy and absolute of position + insuficientAccuracy = false; // will be true if at least one parameter changed too much + for (int iclu = 0; iclu < nMax; iclu++) { + double oldX = xMax[iclu]; + double oldY = yMax[iclu]; + double oldZ = zMax[iclu]; + double oldE = eMax[iclu]; + // new energy, need for weight + eMax[iclu] = 0; + for (int idig = 0; idig < mult; idig++) { + eMax[iclu] += eInClusters[idig][iclu]; + } + xMax[iclu] = 0; + yMax[iclu] = 0; + zMax[iclu] = 0; + double wtot = 0.; + for (int idig = 0; idig < mult; idig++) { + double w = std::max(std::log(eInClusters[idig][iclu] / eMax[iclu]) + mLogWeight, 0.); + xMax[iclu] += x[idig] * w; + yMax[iclu] += y[idig] * w; + zMax[iclu] += z[idig] * w; + wtot += w; + } + if (wtot > 0.) { + xMax[iclu] /= wtot; + yMax[iclu] /= wtot; + zMax[iclu] /= wtot; + } + // Compare variation of parameters + insuficientAccuracy += (std::abs(eMax[iclu] - oldE) > mUnfogingEAccuracy); + insuficientAccuracy += (std::abs(xMax[iclu] - oldX) > mUnfogingXZAccuracy); + insuficientAccuracy += (std::abs(yMax[iclu] - oldY) > mUnfogingXZAccuracy); + insuficientAccuracy += (std::abs(zMax[iclu] - oldZ) > mUnfogingXZAccuracy); + } + nIterations++; + } + + // Iterations finished, add new clusters + for (int iclu = 0; iclu < nMax; iclu++) { + auto& clu = unfoldedClusters.emplace_back(); + clu.setNLM(nMax); + for (int idig = 0; idig < mult; idig++) { + int jdigit = iniClu->getDigitIndex(idig); + int towerId = iniClu->getDigitTowerId(idig); + clu.addDigit(jdigit, towerId, eInClusters[idig][iclu]); + } + } +} + +//============================================================================== +void Clusterizer::evalClusters(std::vector& clusters) +{ + auto& geo = Geometry::instance(); + for (auto& cluster : clusters) { + double x = 0; + double y = 0; + double z = 0; + double wtot = 0; + double etot = cluster.getE(); + for (size_t i = 0; i < cluster.getMultiplicity(); i++) { + float energy = cluster.getDigitEnergy(i); + int towerId = cluster.getDigitTowerId(i); + double xi, yi, zi; + geo.detIdToGlobalPosition(towerId, xi, yi, zi); + double w = std::max(0., mLogWeight + std::log(energy / etot)); + x += w * xi; + y += w * yi; + z += w * zi; + wtot += w; + } + if (wtot != 0) { + x /= wtot; + y /= wtot; + z /= wtot; + } + cluster.setX(x); + cluster.setY(y); + cluster.setZ(z); + + // cluster shape + float chi2 = 0; + int ndf = 0; + float ee = cluster.getE(); + for (size_t i = 0; i < cluster.getMultiplicity(); i++) { + float energy = cluster.getDigitEnergy(i); + int towerId = cluster.getDigitTowerId(i); + double xi, yi, zi; + geo.detIdToGlobalPosition(towerId, xi, yi, zi); + double r = std::sqrt((x - xi) * (x - xi) + (y - yi) * (y - yi) + (z - zi) * (z - zi)); + if (r > 2.2) + continue; + double frac = fCrystalShowerShape->Eval(r); + double rms = fCrystalRMS->Eval(r); + chi2 += std::pow((energy / ee - frac) / rms, 2.); + ndf++; + } + cluster.setChi2(chi2 / ndf); + + // correct cluster energy and z position + float eta = std::abs(cluster.getEta()); + float eCor = 1; + float zCor = 0; + bool isCrystal = geo.isCrystal(cluster.getDigitTowerId(0)); + if (isCrystal) { + eCor = 0.00444 * std::pow(ee, -1.322) + (1.021 + 0.0018 * eta); + if (mApplyCorrectionE) + ee *= eCor; + if (mApplyCorrectionZ) + zCor = (-0.00518682 + 0.730052 * eta - 0.73817 * eta * eta); + } else { + eCor = 0.0033 * std::pow(ee, -2.09) + (1.007 + 0.0667 * eta - 0.108 * eta * eta + 0.0566 * eta * eta * eta); + if (mApplyCorrectionE) + ee *= eCor; + if (mApplyCorrectionZ) + zCor = (-2.13679 + 6.40009 * eta - 3.34233 * eta * eta) + (-0.136425 + 0.401887 * eta - 0.196851 * eta * eta) * ee + (0.00822276 - 0.0242512 * eta + 0.0118986 * eta * eta) * ee * ee; + } + + cluster.setE(ee); + cluster.setZ(cluster.getZ() - zCor); + + // check if cluster is at the edge of detector module + bool isEdge = 0; + for (size_t i = 0; i < cluster.getMultiplicity(); i++) { + int towerId = cluster.getDigitTowerId(i); + if (!geo.isAtTheEdge(towerId)) + continue; + isEdge = 1; + break; + } + cluster.setEdgeFlag(isEdge); + + LOGF(debug, "Cluster coordinates: (%6.2f,%6.2f,%6.2f), eCor=%6.2f zCor=%6.2f", cluster.getX(), cluster.getY(), cluster.getZ(), eCor, zCor); + } +} + +//============================================================================== +int Clusterizer::getNumberOfLocalMax(Cluster& clu, int* maxAt, float* maxAtEnergy) +{ + // Based on MpdEmcClusterizerKI::GetNumberOfLocalMax by D. Peresunko + // Calculates the number of local maxima in the cluster using LocalMaxCut as the minimum + // energy difference between maximum and surrounding digits + auto& geo = Geometry::instance(); + int n = clu.getMultiplicity(); + bool isCrystal = geo.isCrystal(clu.getDigitTowerId(0)); + bool* isLocalMax = new bool[n]; + + for (int i = 0; i < n; i++) { + isLocalMax[i] = false; + float en1 = clu.getDigitEnergy(i); + if (en1 > mClusteringThreshold) + isLocalMax[i] = true; + } + + for (int i = 0; i < n; i++) { + int detId1 = clu.getDigitTowerId(i); + float en1 = clu.getDigitEnergy(i); + for (int j = i + 1; j < n; j++) { + int detId2 = clu.getDigitTowerId(j); + float en2 = clu.getDigitEnergy(j); + if (geo.areNeighboursVertex(detId1, detId2) == 1) { + if (en1 > en2) { + isLocalMax[j] = false; + // but may be digit too is not local max ? + if (en2 > en1 - mLocalMaximumCut) { + isLocalMax[i] = false; + } + } else { + isLocalMax[i] = false; + // but may be digitN is not local max too? + if (en1 > en2 - mLocalMaximumCut) { + isLocalMax[j] = false; + } + } + } // if neighbours + } // digit j + } // digit i + + int iDigitN = 0; + for (int i = 0; i < n; i++) { + if (isLocalMax[i]) { + maxAt[iDigitN] = i; + maxAtEnergy[iDigitN] = clu.getDigitEnergy(i); + iDigitN++; + if (iDigitN >= mNLMMax) { // Note that size of output arrays is limited: + LOGP(error, "Too many local maxima, cluster multiplicity {} region={}", n, isCrystal ? "crystal" : "sampling"); + return 0; + } + } + } + delete[] isLocalMax; + return iDigitN; +} + +//============================================================================== +double Clusterizer::showerShape(double dx, double dz, bool isCrystal) +{ + double x = std::sqrt(dx * dx + dz * dz); + return isCrystal ? fCrystalShowerShape->Eval(x) : fSamplingShowerShape->Eval(x); +} diff --git a/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/ECalReconstructionLinkDef.h b/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/ECalReconstructionLinkDef.h new file mode 100644 index 0000000000000..d69cd8164e717 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/reconstruction/src/ECalReconstructionLinkDef.h @@ -0,0 +1,20 @@ +// 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::ecal::Clusterizer + ; + +#endif diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/CMakeLists.txt b/Detectors/Upgrades/ALICE3/ECal/simulation/CMakeLists.txt index 8c8c5a6bba15f..83de48e38db3a 100644 --- a/Detectors/Upgrades/ALICE3/ECal/simulation/CMakeLists.txt +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/CMakeLists.txt @@ -11,8 +11,13 @@ o2_add_library(ECalSimulation SOURCES src/Detector.cxx + src/Digitizer.cxx PUBLIC_LINK_LIBRARIES O2::ECalBase - O2::ITSMFTSimulation) + O2::DataFormatsECal) o2_target_root_dictionary(ECalSimulation - HEADERS include/ECalSimulation/Detector.h) \ No newline at end of file + HEADERS include/ECalSimulation/Detector.h + include/ECalSimulation/Digitizer.h + ) + +o2_data_file(COPY data DESTINATION Detectors/ECL/simulation) \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/data/simcuts.dat b/Detectors/Upgrades/ALICE3/ECal/simulation/data/simcuts.dat new file mode 100644 index 0000000000000..81aa69990f222 --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/data/simcuts.dat @@ -0,0 +1,14 @@ +* ECAL +* ==== +* +* Med GAM ELEC NHAD CHAD MUON EBREM MUHAB EDEL MUDEL MUPA ANNI BREM COMP DCAY DRAY HADR LOSS MULS PAIR PHOT RAYL STRA +* Air +ECL 0 5.e-5 1.e-4 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 -1. 1 1 1 1 1 1 1 1 1 1 1 -1 +* Lead +ECL 1 1.e-4 1.e-4 1.e-3 1.e-3 1.e-3 1.e-4 1.e-4 1.e-4 1.e-4 -1. 1 1 1 1 1 1 3 1 1 1 1 -1 +* Scintillator +ECL 2 1.e-4 1.e-4 1.e-3 1.e-3 1.e-3 1.e-4 1.e-4 1.e-4 1.e-4 -1. 1 1 1 1 1 1 3 1 1 1 1 -1 +* Aluminium +ECL 3 1.e-4 1.e-4 1.e-3 1.e-3 1.e-3 1.e-4 1.e-4 1.e-4 1.e-4 -1. 1 1 1 1 1 1 3 1 1 1 1 -1 +* PWO +ECL 4 5.e-5 1.e-4 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 1.e-5 -1. 1 1 1 1 1 1 1 1 1 1 1 -1 diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Detector.h b/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Detector.h index 14664092a8718..849dd69f85f2b 100644 --- a/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Detector.h +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Detector.h @@ -8,18 +8,18 @@ // 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. -// -// Design and equations: Nicola Nicassio nicola.nicassio@cern.ch -// + +/// \file Detector.h +/// \brief ECal geometry creation and hit processing +/// +/// \author Evgeny Kryshen #ifndef ALICEO2_ECAL_DETECTOR_H #define ALICEO2_ECAL_DETECTOR_H -#include "DetectorsBase/Detector.h" -#include "ITSMFTSimulation/Hit.h" - -#include "ECalBase/GeometryTGeo.h" - +#include +#include +#include #include #include @@ -27,62 +27,35 @@ namespace o2 { namespace ecal { - class Detector : public o2::base::DetImpl { public: - Detector(bool active); - Detector(); + Detector(bool active = 1); ~Detector(); - void ConstructGeometry() override; - - o2::itsmft::Hit* addHit(int trackID, int detID, const TVector3& startPos, const TVector3& endPos, - const TVector3& startMom, double startE, double endTime, double eLoss, - unsigned char startStatus, unsigned char endStatus); - // Mandatory overrides - void BeginPrimary() override { ; } - void FinishPrimary() override { ; } + void ConstructGeometry() override; + void BeginPrimary() override {} + void FinishPrimary() override {} void InitializeO2Detector() override; - void PostTrack() override { ; } - void PreTrack() override { ; } + void PostTrack() override {} + void PreTrack() override {} bool ProcessHits(FairVolume* v = nullptr) override; - void EndOfEvent() override; + void EndOfEvent() override { Reset(); } void Register() override; void Reset() override; + std::vector* getHits(int iColl) const { return !iColl ? mHits : nullptr; } - // Custom memer functions - std::vector* getHits(int iColl) const - { - if (!iColl) { - return mHits; - } - return nullptr; - } - + private: void createMaterials(); void createGeometry(); - - private: - // Transient data about track passing the sensor - struct TrackData { - bool mHitStarted; // hit creation started - unsigned char mTrkStatusStart; // track status flag - TLorentzVector mPositionStart; // position at entrance - TLorentzVector mMomentumStart; // momentum - double mEnergyLoss; // energy loss - } mTrackData; //! transient data - - GeometryTGeo* mGeometryTGeo; //! - std::vector* mHits; // ITSMFT ones for the moment - - void defineSensitiveVolumes(); - float mInnerRadius; - float mOuterRadius; - float mLength; - - bool mEnableEndcap{true}; + void defineSamplingFactor(); + std::unordered_map mSuperParentIndices; //! Super parent indices (track index - superparent index) + int currentTrackId = -1; // current track index + int superparentId = -1; // superparent index + GeometryTGeo* mGeometryTGeo; //! + std::vector* mHits; //! + double mSamplingFactorTransportModel = 1.; protected: template @@ -104,4 +77,5 @@ struct UseShm { } // namespace base } // namespace o2 #endif -#endif \ No newline at end of file + +#endif // ALICEO2_ECAL_DETECTOR_H \ No newline at end of file diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Digitizer.h b/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Digitizer.h new file mode 100644 index 0000000000000..91213fa90b63a --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/include/ECalSimulation/Digitizer.h @@ -0,0 +1,58 @@ +// 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. + +/// \file Digitizer.h +/// \brief Digitization of ECal MC information +/// +/// \author Evgeny Kryshen + +#ifndef ALICEO2_ECAL_DIGITIZER_H +#define ALICEO2_ECAL_DIGITIZER_H +#include +#include +#include +#include +#include + +using o2::ecal::Digit; +using o2::ecal::Hit; +using o2::ecal::MCLabel; + +namespace o2 +{ +namespace ecal +{ +class Digitizer +{ + public: + Digitizer(); + ~Digitizer() = default; + Digitizer(const Digitizer&) = delete; + Digitizer& operator=(const Digitizer&) = delete; + void init() {} + void finish() {} + void processHits(const std::vector* mHits, std::vector& digits, o2::dataformats::MCTruthContainer& labels, int collId); + void setThreshold(double threshold) { mThreshold = threshold; } + void setSmearCrystal(bool smearCrystal) { mSmearCrystal = smearCrystal; } + void setSamplingFraction(double fraction) { mSamplingFraction = fraction; } + void setCrystalPePerGeV(double pePerGeV) { mCrystalPePerGeV = pePerGeV; } + + private: + std::vector mArrayD; + bool mSmearCrystal = 0; + double mThreshold = 0.001; + double mSamplingFraction = 9.8; + double mCrystalPePerGeV = 4000; +}; +} // namespace ecal +} // namespace o2 + +#endif // ALICEO2_ECAL_DIGITIZER_H diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx b/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx index aeb58649fa4c5..93089bb8ced14 100644 --- a/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/src/Detector.cxx @@ -9,45 +9,40 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include +/// \file Detector.cxx +/// \brief ECal geometry creation and hit processing +/// +/// \author Evgeny Kryshen +#include #include #include +#include #include #include -#include - -#include "DetectorsBase/Stack.h" -#include "ITSMFTSimulation/Hit.h" -#include "ECalSimulation/Detector.h" -#include "ECalBase/ECalBaseParam.h" - -using o2::itsmft::Hit; +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using o2::ecal::Hit; namespace o2 { namespace ecal { - -Detector::Detector() - : o2::base::DetImpl("ECL", true), - mTrackData(), - mHits(o2::utils::createSimVector()) -{ -} - Detector::Detector(bool active) - : o2::base::DetImpl("ECL", true), - mTrackData(), - mHits(o2::utils::createSimVector()) + : o2::base::DetImpl("ECL", active), + mHits(o2::utils::createSimVector()) { - auto& ecalPars = ECalBaseParam::Instance(); - mInnerRadius = ecalPars.rMin; - mOuterRadius = ecalPars.rMax; - mLength = ecalPars.zLength; - mEnableEndcap = ecalPars.enableFwdEndcap; } +//============================================================================== Detector::~Detector() { if (mHits) { @@ -55,190 +50,344 @@ Detector::~Detector() } } +//============================================================================== void Detector::ConstructGeometry() { createMaterials(); createGeometry(); + defineSamplingFactor(); } +//============================================================================== +void Detector::defineSamplingFactor() +{ + TString mcname = TVirtualMC::GetMC()->GetName(); + TString mctitle = TVirtualMC::GetMC()->GetTitle(); + LOGP(info, "Defining sampling factor for mc={}' and title='{}'", mcname.Data(), mctitle.Data()); + if (mcname.Contains("Geant3")) { + mSamplingFactorTransportModel = 0.983; + } else if (mcname.Contains("Geant4")) { + mSamplingFactorTransportModel = 1.; + } +} + +//============================================================================== void Detector::createMaterials() { - int ifield = 2; // ? - float fieldm = 10.0; // ? + LOGP(info, "Creating materials for ECL"); + + // Air + float aAir[4] = {12.0107, 14.0067, 15.9994, 39.948}; + float zAir[4] = {6., 7., 8., 18.}; + float wAir[4] = {0.000124, 0.755267, 0.231781, 0.012827}; + float dAir = 1.20479E-3; + Mixture(0, "Air", aAir, zAir, dAir, 4, wAir); + + // Pb + Material(1, "Pb", 207.2, 82, 11.35, 0.56, 0., nullptr, 0); + + // Polysterene scintillator (CH) + float aP[2] = {12.011, 1.00794}; + float zP[2] = {6.0, 1.0}; + float wP[2] = {1.0, 1.0}; + float dP = 1.032; + Mixture(2, "Scintillator", aP, zP, dP, -2, wP); + + // Al + Material(3, "Al", 26.98, 13., 2.7, 8.9, 999., nullptr, 0); + + // PWO crystals + float aX[3] = {207.19, 183.85, 16.0}; + float zX[3] = {82.0, 74.0, 8.0}; + float wX[3] = {1.0, 1.0, 4.0}; + float dX = 8.28; + Mixture(4, "PbWO4", aX, zX, dX, -3, wX); + + int ifield = 2; // magnetic field flag + float fieldm = 10.; // maximum field value (in Kilogauss) + float deemax = 0.1; // maximum fractional energy loss in one step (0 < deemax <=1) + float tmaxfd = 10.; // maximum angle due to field permitted in one step (in degrees) o2::base::Detector::initFieldTrackingParams(ifield, fieldm); - - float tmaxfdLead = 0.1; // .10000E+01; // Degree - float stemaxLead = .10000E+01; // cm - float deemaxLead = 0.1; // 0.30000E-02; // Fraction of particle's energy 0clear(); + } + mSuperParentIndices.clear(); + currentTrackId = -1; + superparentId = -1; } -void Detector::EndOfEvent() { Reset(); } - +//============================================================================== void Detector::Register() { // This will create a branch in the output tree called Hit, setting the last // parameter to kFALSE means that this collection will not be written to the file, // it will exist only during the simulation - + LOGP(info, "Registering hits"); if (FairRootManager::Instance()) { FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, true); } } +//============================================================================== void Detector::createGeometry() { LOGP(info, "Creating ECal geometry"); - TGeoManager* geoManager = gGeoManager; - TGeoVolume* vALIC = geoManager->GetVolume("barrel"); + TGeoVolume* vALIC = gGeoManager->GetVolume("barrel"); if (!vALIC) { LOGP(fatal, "Could not find barrel volume while constructing ECal geometry"); } new TGeoVolumeAssembly(GeometryTGeo::getECalVolPattern()); - TGeoVolume* vECal = geoManager->GetVolume(GeometryTGeo::getECalVolPattern()); + TGeoVolume* vECal = gGeoManager->GetVolume(GeometryTGeo::getECalVolPattern()); vALIC->AddNode(vECal, 2, new TGeoTranslation(0, 30., 0)); + vECal->SetTitle("ECalVol"); + + TGeoMedium* medAir = gGeoManager->GetMedium("ECL_Air"); + TGeoMedium* medPb = gGeoManager->GetMedium("ECL_Pb"); + TGeoMedium* medAl = gGeoManager->GetMedium("ECL_Al"); + TGeoMedium* medSc = gGeoManager->GetMedium("ECL_Scintillator"); + TGeoMedium* medPWO = gGeoManager->GetMedium("ECL_Crystal"); + + // Get relevant parameters + auto& pars = ECalBaseParam::Instance(); + auto& geo = Geometry::instance(); + + double rMin = pars.rMin; + double rMax = pars.rMax; + double layerThickness = pars.pbLayerThickness + pars.scLayerThickness; + double samplingModL = pars.frontPlateThickness + layerThickness * pars.nSamplingLayers - pars.pbLayerThickness; + double crystalAlpha = geo.getCrystalAlpha(); + double samplingAlpha = geo.getSamplingAlpha(); + double tanCrystalAlpha = std::tan(crystalAlpha); + double tanSamplingAlpha = std::tan(samplingAlpha); + + double sectorL = rMax - rMin; + double crystalThetaMin = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ - 1) - crystalAlpha; + double crystalHlMin = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ - 1); + double crystalHlMax = crystalHlMin + sectorL / std::tan(crystalThetaMin); + double crystalHwMin = geo.getCrystalModW() / 2.; + double crystalHwMax = crystalHwMin * rMax / rMin; + auto crystalSectorShape = new TGeoTrap(sectorL / 2., 0, 0, crystalHlMin, crystalHwMin, crystalHwMin, 0, crystalHlMax, crystalHwMax, crystalHwMax, 0); + auto crystalSectorVolume = new TGeoVolume("crystalSectorVolume", crystalSectorShape, medAir); + AddSensitiveVolume(crystalSectorVolume); + crystalSectorVolume->SetLineColor(kCyan + 1); + crystalSectorVolume->SetTransparency(0); + + double samplingThetaAtMinZ = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ) + samplingAlpha; + double samplingThetaAtMaxZ = geo.getFrontFaceCenterTheta(pars.nCrystalModulesZ + pars.nSamplingModulesZ - 1) - samplingAlpha; + double samplingMinZatMinR = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ) - geo.getSamplingModW() / std::sin(samplingThetaAtMinZ); + double samplingMaxZatMinR = geo.getFrontFaceZatMinR(pars.nCrystalModulesZ + pars.nSamplingModulesZ - 1); + double samplingMinZatMaxR = samplingMinZatMinR + sectorL / std::tan(samplingThetaAtMinZ); + double samplingMaxZatMaxR = samplingMaxZatMinR + sectorL / std::tan(samplingThetaAtMaxZ); + double hlMin = (samplingMaxZatMinR - samplingMinZatMinR) / 2.; + double hlMax = (samplingMaxZatMaxR - samplingMinZatMaxR) / 2.; + double zCenterMin = (samplingMaxZatMinR + samplingMinZatMinR) / 2.; + double zCenterMax = (samplingMaxZatMaxR + samplingMinZatMaxR) / 2.; + double zCenter = (zCenterMax + zCenterMin) / 2.; + double thetaCenter = std::atan((zCenterMax - zCenterMin) / sectorL) * TMath::RadToDeg(); + double samplingHwMin = geo.getSamplingModW() / 2.; + double samplingHwMax = samplingHwMin * rMax / rMin; + auto samplingSectorShape = new TGeoTrap(sectorL / 2., thetaCenter, 90, hlMin, samplingHwMin, samplingHwMin, 0, hlMax, samplingHwMax, samplingHwMax, 0); + auto samplingSectorVolume = new TGeoVolume("samplingSectorVolume", samplingSectorShape, medAir); + AddSensitiveVolume(samplingSectorVolume); + samplingSectorVolume->SetLineColor(kBlue + 1); + samplingSectorVolume->SetTransparency(0); + + double sectorR = rMin + sectorL / 2.; + for (int ism = 0; ism < pars.nSuperModules; ism++) { + // crystal + for (int i = 0; i < pars.nCrystalModulesPhi; i++) { + int row = ism * pars.nCrystalModulesPhi + i; + double phi0 = geo.getFrontFaceCenterCrystalPhi(row); + double x = sectorR * std::cos(phi0); + double y = sectorR * std::sin(phi0); + auto rot = new TGeoRotation(Form("ecalcrystalsecrot%d", row), 90 + phi0 * TMath::RadToDeg(), 90, 0); + vECal->AddNode(crystalSectorVolume, row, new TGeoCombiTrans(x, y, 0., rot)); + } + // sampling + for (int i = 0; i < pars.nSamplingModulesPhi; i++) { + int row = ism * pars.nSamplingModulesPhi + i; + double phi0 = geo.getFrontFaceCenterSamplingPhi(row); + double x = sectorR * std::cos(phi0); + double y = sectorR * std::sin(phi0); + auto rot1 = new TGeoRotation(Form("ecalsamplingsec1rot%d", row), 90 + phi0 * TMath::RadToDeg(), 90, 0.); + auto rot2 = new TGeoRotation(Form("ecalsamplingsec2rot%d", row), 90 + phi0 * TMath::RadToDeg(), 90, 180); + vECal->AddNode(samplingSectorVolume, 2 * row + 0, new TGeoCombiTrans(x, y, zCenter, rot1)); + vECal->AddNode(samplingSectorVolume, 2 * row + 1, new TGeoCombiTrans(x, y, -zCenter, rot2)); + } + } - char vstrng[100] = "ECalVol"; - vECal->SetTitle(vstrng); + for (int m = 0; m < pars.nCrystalModulesZ; m++) { + double tanBeta = geo.getTanBeta(m); + double dx1 = crystalHwMin; + double dx2 = crystalHwMin + pars.crystalModuleLength * tanCrystalAlpha; + double dy1 = crystalHwMin; + double dy2 = crystalHwMin + pars.crystalModuleLength * tanBeta; + double dz = pars.crystalModuleLength / 2.; + auto crystalModuleShape = new TGeoTrd2(dx1, dx2, dy1, dy2, dz); + auto crystalModuleVolume = new TGeoVolume(Form("crystalmodule%d", m), crystalModuleShape, medPWO); + AddSensitiveVolume(crystalModuleVolume); + crystalModuleVolume->SetLineColor(kCyan + 1); + crystalModuleVolume->SetTransparency(0); + double theta = geo.getFrontFaceCenterTheta(m); + double r = geo.getFrontFaceCenterR(m); + double z = geo.getFrontFaceCenterZ(m); + ROOT::Math::XYPoint pFrontFace(z, r - sectorR); + ROOT::Math::Polar2DVector vFrontFaceToCenter(dz, theta); + ROOT::Math::XYPoint pc = pFrontFace + vFrontFaceToCenter; + auto rot1 = new TGeoRotation(Form("ecalcrystalrot%d", 2 * m), 0, 270 + theta * TMath::RadToDeg(), 90); + crystalSectorVolume->AddNode(crystalModuleVolume, 2 * m, new TGeoCombiTrans(0, pc.x(), pc.y(), rot1)); + auto rot2 = new TGeoRotation(Form("ecalcrystalrot%d", 2 * m + 1), 0, 90 - theta * TMath::RadToDeg(), 90); + crystalSectorVolume->AddNode(crystalModuleVolume, 2 * m + 1, new TGeoCombiTrans(0, -pc.x(), pc.y(), rot2)); + } - // Build the ECal cylinder - auto& matmgr = o2::base::MaterialManager::Instance(); - TGeoMedium* medPb = matmgr.getTGeoMedium("ECL_LEAD"); - TGeoTube* ecalShape = new TGeoTube("ECLsh", mInnerRadius, mOuterRadius, mLength); - TGeoVolume* ecalVol = new TGeoVolume("ECL", ecalShape, medPb); - ecalVol->SetLineColor(kAzure - 9); - ecalVol->SetTransparency(0); - vECal->AddNode(ecalVol, 1, nullptr); + for (int m = 0; m < pars.nSamplingModulesZ; m++) { + int k = pars.nCrystalModulesZ + m; + double tanBeta = geo.getTanBeta(k); + double dx1 = samplingHwMin; + double dx2 = samplingHwMin + samplingModL * tanSamplingAlpha; + double dy1 = samplingHwMin; + double dy2 = samplingHwMin + samplingModL * tanBeta; + double dz = samplingModL / 2.; + auto samplingModuleShape = new TGeoTrd2(dx1, dx2, dy1, dy2, dz); + auto samplingModuleVolume = new TGeoVolume(Form("samplingmodule%d", m), samplingModuleShape, medSc); + AddSensitiveVolume(samplingModuleVolume); + samplingModuleVolume->SetLineColor(kAzure - 9); + samplingModuleVolume->SetTransparency(0); + double theta = geo.getFrontFaceCenterTheta(k); + double r = geo.getFrontFaceCenterR(k); + double z = geo.getFrontFaceCenterZ(k); + ROOT::Math::XYPoint pFrontFace(z - zCenter, r - sectorR); + ROOT::Math::Polar2DVector vFrontFaceToCenter(dz, theta); + ROOT::Math::XYPoint pc = pFrontFace + vFrontFaceToCenter; + auto rot1 = new TGeoRotation(Form("ecalsamplingrot%d", m), 0, 270 + theta * TMath::RadToDeg(), 90); + samplingSectorVolume->AddNode(samplingModuleVolume, m, new TGeoCombiTrans(0, pc.x(), pc.y(), rot1)); + + // adding front aluminium plate into the volume + double fdx1 = dx1; + double fdx2 = dx1 + pars.frontPlateThickness * tanSamplingAlpha; + double fdy1 = dy1; + double fdy2 = fdy1 + pars.frontPlateThickness * tanBeta; + double fdz = pars.frontPlateThickness / 2.; + auto frontShape = new TGeoTrd2(fdx1, fdx2, fdy1, fdy2, fdz); + auto frontVolume = new TGeoVolume(Form("front%d", m), frontShape, medAl); + samplingModuleVolume->AddNode(frontVolume, 0, new TGeoTranslation(0., 0., -dz + pars.frontPlateThickness / 2.)); + AddSensitiveVolume(frontVolume); + frontVolume->SetLineColor(kAzure - 7); + frontVolume->SetTransparency(0); + + // adding lead plates + for (int i = 0; i < pars.nSamplingLayers - 1; i++) { + double lz1 = pars.frontPlateThickness + pars.scLayerThickness + layerThickness * i; + double lz2 = lz1 + pars.pbLayerThickness; + double lzc = -dz + (lz1 + lz2) / 2.; + double ldx1 = dx1 + lz1 * tanSamplingAlpha; + double ldx2 = dx1 + lz2 * tanSamplingAlpha; + double ldy1 = dy1 + lz1 * tanBeta; + double ldy2 = dy1 + lz2 * tanBeta; + double ldz = pars.pbLayerThickness / 2.; + auto leadShape = new TGeoTrd2(ldx1, ldx2, ldy1, ldy2, ldz); + auto leadVolume = new TGeoVolume(Form("lead%d_%d", m, i), leadShape, medPb); + samplingModuleVolume->AddNode(leadVolume, i, new TGeoTranslation(0., 0., lzc)); + AddSensitiveVolume(leadVolume); + leadVolume->SetLineColor(kAzure - 7); + leadVolume->SetTransparency(0); + } + } - if (mEnableEndcap) { + if (pars.enableFwdEndcap) { // Build the ecal endcap - TGeoTube* ecalEndcapShape = new TGeoTube("ECLECsh", 15.f, 160.f, 0.5 * (mOuterRadius - mInnerRadius)); + TGeoTube* ecalEndcapShape = new TGeoTube("ECLECsh", 15.f, 160.f, 0.5 * (rMax - rMin)); TGeoVolume* ecalEndcapVol = new TGeoVolume("ECLEC", ecalEndcapShape, medPb); ecalEndcapVol->SetLineColor(kAzure - 9); ecalEndcapVol->SetTransparency(0); vECal->AddNode(ecalEndcapVol, 1, new TGeoTranslation(0, 0, -450.f)); } + // gGeoManager->CloseGeometry(); + // gGeoManager->CheckOverlaps(0.0001); } -void Detector::Reset() -{ - if (!o2::utils::ShmManager::Instance().isOperational()) { - mHits->clear(); - } -} - +//============================================================================== bool Detector::ProcessHits(FairVolume* vol) { - // This method is called from the MC stepping - if (!(fMC->TrackCharge())) { - return false; - } - - int lay = vol->getVolumeId(); - int volID = vol->getMCid(); - - // Is it needed to keep a track reference when the outer ITS volume is encountered? + LOGP(debug, "Processing hits"); auto stack = (o2::data::Stack*)fMC->GetStack(); - if (fMC->IsTrackExiting() && (lay == 0)) { - o2::TrackReference tr(*fMC, GetDetId()); - tr.setTrackID(stack->GetCurrentTrackNumber()); - tr.setUserId(lay); - stack->addTrackReference(tr); - } - bool startHit = false, stopHit = false; - unsigned char status = 0; - if (fMC->IsTrackEntering()) { - status |= Hit::kTrackEntering; - } - if (fMC->IsTrackInside()) { - status |= Hit::kTrackInside; - } - if (fMC->IsTrackExiting()) { - status |= Hit::kTrackExiting; - } - if (fMC->IsTrackOut()) { - status |= Hit::kTrackOut; - } - if (fMC->IsTrackStop()) { - status |= Hit::kTrackStopped; - } - if (fMC->IsTrackAlive()) { - status |= Hit::kTrackAlive; + int trackId = stack->GetCurrentTrackNumber(); + int parentId = stack->GetCurrentParentTrackNumber(); + + if (trackId != currentTrackId) { + auto superparentIndexIt = mSuperParentIndices.find(parentId); + if (superparentIndexIt != mSuperParentIndices.end()) { + superparentId = superparentIndexIt->second; + mSuperParentIndices[trackId] = superparentIndexIt->second; + } else { + // for new incoming tracks the superparent index is equal to the track ID (for recursion) + mSuperParentIndices[trackId] = trackId; + superparentId = trackId; + } + currentTrackId = trackId; } - // track is entering or created in the volume - if ((status & Hit::kTrackEntering) || (status & Hit::kTrackInside && !mTrackData.mHitStarted)) { - startHit = true; - } else if ((status & (Hit::kTrackExiting | Hit::kTrackOut | Hit::kTrackStopped))) { - stopHit = true; + double eloss = fMC->Edep(); + if (eloss < DBL_EPSILON) { + return false; } - // increment energy loss at all steps except entrance - if (!startHit) { - mTrackData.mEnergyLoss += fMC->Edep(); - } - if (!(startHit | stopHit)) { - return false; // do noting + TString volName = vol->GetName(); + bool isCrystal = volName.Contains("crystalmodule"); + bool isSampling = volName.Contains("samplingmodule"); + + if (!isCrystal && !isSampling) { + return false; } - if (startHit) { - mTrackData.mEnergyLoss = 0.; - fMC->TrackMomentum(mTrackData.mMomentumStart); - fMC->TrackPosition(mTrackData.mPositionStart); - mTrackData.mTrkStatusStart = status; - mTrackData.mHitStarted = true; + if (isCrystal) + LOGP(debug, "Processing crystal {}", volName.Data()); + else { + eloss *= mSamplingFactorTransportModel; + LOGP(debug, "Processing scintillator {}", volName.Data()); } - if (stopHit) { - TLorentzVector positionStop; - fMC->TrackPosition(positionStop); - // Retrieve the indices with the volume path - int stave(0), halfstave(0), chipinmodule(0), module; - fMC->CurrentVolOffID(1, chipinmodule); - fMC->CurrentVolOffID(2, module); - fMC->CurrentVolOffID(3, halfstave); - fMC->CurrentVolOffID(4, stave); - - Hit* p = addHit(stack->GetCurrentTrackNumber(), lay, mTrackData.mPositionStart.Vect(), positionStop.Vect(), - mTrackData.mMomentumStart.Vect(), mTrackData.mMomentumStart.E(), positionStop.T(), - mTrackData.mEnergyLoss, mTrackData.mTrkStatusStart, status); - // p->SetTotalEnergy(vmc->Etot()); - - // RS: not sure this is needed - // Increment number of Detector det points in TParticle + int sectorId, moduleId; + fMC->CurrentVolID(moduleId); + fMC->CurrentVolOffID(1, sectorId); + int cellID = Geometry::instance().getCellID(moduleId, sectorId, isCrystal); + LOGP(debug, "isCrystal={} sectorId={} moduleId={} cellID={} eloss={}", isCrystal, sectorId, moduleId, cellID, eloss); + + int trackID = superparentId; + auto hit = std::find_if(mHits->begin(), mHits->end(), [cellID, trackID](const Hit& hit) { return hit.GetTrackID() == trackID && hit.GetCellID() == cellID; }); + if (hit == mHits->end()) { + float posX, posY, posZ, momX, momY, momZ, energy; + fMC->TrackPosition(posX, posY, posZ); + fMC->TrackMomentum(momX, momY, momZ, energy); + auto pos = math_utils::Point3D(posX, posY, posZ); + auto mom = math_utils::Vector3D(momX, momY, momZ); + float time = fMC->TrackTime() * 1e9; // time in ns + mHits->emplace_back(trackID, cellID, pos, mom, time, eloss); stack->addHit(GetDetId()); + } else { + hit->SetEnergyLoss(hit->GetEnergyLoss() + eloss); } - return true; } -o2::itsmft::Hit* Detector::addHit(int trackID, int detID, const TVector3& startPos, const TVector3& endPos, - const TVector3& startMom, double startE, double endTime, double eLoss, unsigned char startStatus, - unsigned char endStatus) -{ - mHits->emplace_back(trackID, detID, startPos, endPos, startMom, startE, endTime, eLoss, startStatus, endStatus); - return &(mHits->back()); -} - } // namespace ecal } // namespace o2 -ClassImp(o2::ecal::Detector); \ No newline at end of file +ClassImp(o2::ecal::Detector); diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/src/Digitizer.cxx b/Detectors/Upgrades/ALICE3/ECal/simulation/src/Digitizer.cxx new file mode 100644 index 0000000000000..f213ba563d86d --- /dev/null +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/src/Digitizer.cxx @@ -0,0 +1,89 @@ +// 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. + +/// \file Digitizer.cxx +/// \brief Digitization of ECal MC information +/// +/// \author Evgeny Kryshen + +#include +#include +#include + +#include +#include +#include +#include + +using namespace o2::ecal; + +//============================================================================== +Digitizer::Digitizer() +{ + auto& geo = Geometry::instance(); + mArrayD.resize(geo.getNrows() * geo.getNcols()); +} + +//============================================================================== +void Digitizer::processHits(const std::vector* hits, std::vector& digits, o2::dataformats::MCTruthContainer& labels, int collId) +{ + digits.clear(); + labels.clear(); + + LOGP(debug, "nHits = {}", hits->size()); + auto& geo = Geometry::instance(); + + for (int i = 0; i < mArrayD.size(); i++) { + mArrayD[i].setAmplitude(0); + mArrayD[i].setTimeStamp(1000); + mArrayD[i].setTower(i); + mArrayD[i].setLabel(-1); + // TODO: simulate noise + } + + for (auto& hit : *hits) { + int cellID = hit.GetCellID(); + double eloss = hit.GetEnergyLoss(); + double t = hit.GetTime(); + double elossSmeared = eloss; + bool isCrystal = geo.isCrystal(cellID); + if (isCrystal) { // crystal + double elossSmearedNpe = gRandom->Poisson(eloss * mCrystalPePerGeV) / mCrystalPePerGeV; + if (mSmearCrystal) + elossSmeared = elossSmearedNpe * gRandom->Gaus(1, 0.007); // light attenuation in crystals + } else { // sampling + elossSmeared *= mSamplingFraction; + } + + Digit& digit = mArrayD[cellID]; + digit.setAmplitude(digit.getAmplitude() + elossSmeared); + if (t < digit.getTimeStamp()) + digit.setTimeStamp(t); // setting earliest time, TODO: add time smearing + LOGF(debug, " crystal: %d cellID = %5d, eloss = %8.5f elossSmeared = %8.5f time = %8.5f", isCrystal, cellID, eloss, elossSmeared, t); + + // Adding MC info + MCLabel label(hit.GetTrackID(), collId, 0, false, hit.GetEnergyLoss()); + int labelIndex = digit.getLabel(); + if (labelIndex == -1) { + labelIndex = labels.getIndexedSize(); + labels.addElement(labelIndex, label); + digit.setLabel(labelIndex); + } else { + labels.addElementRandomAccess(labelIndex, label); + } + } // hits + + for (int i = 0; i < mArrayD.size(); i++) { + if (mArrayD[i].getAmplitude() > mThreshold) { + digits.push_back(mArrayD[i]); + } + } +} diff --git a/Detectors/Upgrades/ALICE3/ECal/simulation/src/ECalSimulationLinkDef.h b/Detectors/Upgrades/ALICE3/ECal/simulation/src/ECalSimulationLinkDef.h index 167342773f196..5d7383f086362 100644 --- a/Detectors/Upgrades/ALICE3/ECal/simulation/src/ECalSimulationLinkDef.h +++ b/Detectors/Upgrades/ALICE3/ECal/simulation/src/ECalSimulationLinkDef.h @@ -17,5 +17,6 @@ #pragma link C++ class o2::ecal::Detector + ; #pragma link C++ class o2::base::DetImpl < o2::ecal::Detector> + ; +#pragma link C++ class o2::ecal::Digitizer + ; #endif