diff --git a/Examples/options/LCIO_read_pan.py b/Examples/options/LCIO_read_pan.py
index 32292b36c7f0263841b7dc7663ebaea681b2dabf..c85a3cae511ceeb1cddb600fb5d05e5af0e5e484 100644
--- a/Examples/options/LCIO_read_pan.py
+++ b/Examples/options/LCIO_read_pan.py
@@ -52,6 +52,9 @@ gearSvc.GearXMLFile = "../Detector/DetCEPCv4/compact/FullDetGear.xml"
from Configurables import PandoraPFAlg
pandoralg = PandoraPFAlg("PandoraPFAlg")
+pandoralg.use_dd4hep_geo = False
+pandoralg.use_dd4hep_decoder = False
+pandoralg.use_preshower = False
pandoralg.collections = [
"MCParticle:MCParticle",
"CalorimeterHit:ECALBarrel",
diff --git a/Examples/options/tut_detsim_pandora.py b/Examples/options/tut_detsim_pandora.py
index af80c1e6565ec12dc9cef76faa1bfa7b8fd62fd8..98fd9617dc55515a3e61249b520367d8354ceaf1 100644
--- a/Examples/options/tut_detsim_pandora.py
+++ b/Examples/options/tut_detsim_pandora.py
@@ -114,11 +114,14 @@ example_CaloDigiAlg.CaloAssociationCollection = "RecoCaloAssociation_ECALBarr
##############################################################################
from Configurables import GearSvc
gearSvc = GearSvc("GearSvc")
-gearSvc.GearXMLFile = "../Detector/DetCEPCv4/compact/FullDetGear.xml"
+gearSvc.GearXMLFile = "Detector/DetCEPCv4/compact/FullDetGear.xml"
##############################################################################
from Configurables import PandoraPFAlg
pandoralg = PandoraPFAlg("PandoraPFAlg")
+pandoralg.use_dd4hep_geo = True
+pandoralg.use_dd4hep_decoder = True
+pandoralg.use_preshower = False
pandoralg.collections = [
"MCParticle:MCParticle",
"CalorimeterHit:ECALBarrel",
@@ -143,7 +146,7 @@ pandoralg.WriteReconstructedParticleCollection = "PandoraPFOs"
pandoralg.WriteVertexCollection = "PandoraPFANewStartVertices"
pandoralg.AnaOutput = "Ana.root"
-pandoralg.PandoraSettingsDefault_xml = "../Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml"
+pandoralg.PandoraSettingsDefault_xml = "Reconstruction/PFA/Pandora/PandoraSettingsDefault.xml"
#### Do not chage the collection name, only add or remove ###############
pandoralg.TrackCollections = ["MarlinTrkTracks"]
pandoralg.ECalCaloHitCollections= ["ECALBarrel", "ECALEndcap", "ECALOther"]
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt b/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt
index be0c48674b7465d7ff1bfd607eff6b172439933c..3842a7aaf8e49e03f6a91de2a98ff5c24b66f71b 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/CMakeLists.txt
@@ -1,10 +1,10 @@
gaudi_subdir(GaudiPandora v0r0)
find_package(LCIO REQUIRED )
+find_package(DD4hep COMPONENTS DDG4 REQUIRED)
find_package(GEAR REQUIRED)
message("ENV GEAR: $ENV{GEAR}")
-find_package(CLHEP REQUIRED;CONFIG)
find_package(EDM4HEP REQUIRED )
include_directories(${EDM4HEP_INCLUDE_DIR})
@@ -16,14 +16,14 @@ include_directories(${LCContent_INCLUDE_DIRS})
link_libraries(${LCContent_LIBRARIES})
-list(APPEND CMAKE_MODULE_PATH "$ENV{ROOTSYS}/etc/cmake/")
-find_package(ROOT 5.26.00 REQUIRED COMPONENTS Eve Geom RGL EG)
+find_package(ROOT COMPONENTS MathCore Physics GenVector Geom REQUIRED)
gaudi_depends_on_subdirs(
Service/EventSeeder
Service/GearSvc
Utilities/DataHelper
+ Detector/DetInterface
)
set(dir_srcs
@@ -33,12 +33,13 @@ set(dir_srcs
src/CaloHitCreator.cpp
src/TrackCreator.cpp
src/PfoCreator.cpp
+ src/Utility.cpp
)
set(dir_include include)
# Modules
gaudi_add_module(GaudiPandora ${dir_srcs}
- INCLUDE_DIRS ${dir_include} GaudiKernel FWCore ${CLHEP_INCLUDE_DIR} ${LCIO_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS} gear
- LINK_LIBRARIES GaudiKernel FWCore ${CLHEP_LIBRARIES} ROOT ${LCIO_LIBRARIES} ${GEAR_LIBRARIES} DataHelperLib
+ INCLUDE_DIRS ${dir_include} GaudiKernel FWCore CLHEP ${LCIO_INCLUDE_DIRS} ROOT gear DD4hep
+ LINK_LIBRARIES GaudiKernel FWCore CLHEP ROOT ${LCIO_LIBRARIES} ${GEAR_LIBRARIES} DataHelperLib DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec
-Wl,--no-as-needed
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
-Wl,--as-needed
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
index 91fd1ccd04f1de4ce7de89d83048171882564ab2..b947f49d25bf7e6147302d4c167e3bac7b026227 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/CaloHitCreator.h
@@ -11,6 +11,14 @@
#include "GaudiKernel/ISvcLocator.h"
#include "edm4hep/CalorimeterHit.h"
+#include "DetInterface/IGeomSvc.h"
+#include "gear/LayerLayout.h"
+#include <DD4hep/DD4hepUnits.h>
+#include <DD4hep/DetType.h>
+#include <DD4hep/DetectorSelector.h>
+#include <DD4hep/Detector.h>
+
+#include "Utility.h"
#include "gear/LayerLayout.h"
#include "Api/PandoraApi.h"
@@ -89,6 +97,9 @@ public:
float m_eCalScToHadGeVBarrel; ///< The calibration from deposited Sc-layer energy on the endcaps to hadronic energy
float m_eCalSiToHadGeVEndCap; ///< The calibration from deposited Si-layer energy on the enecaps to hadronic energy
float m_eCalScToHadGeVEndCap; ///< The calibration from deposited Sc-layer energy on the endcaps to hadronic energy
+ bool m_use_dd4hep_geo; ///
+ bool m_use_dd4hep_decoder; ///
+ bool m_use_preshower; ///
};
/**
@@ -165,6 +176,8 @@ private:
*/
void GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const;
+ void GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const;
/**
* @brief Get barrel specific calo hit properties: cell size, absorber radiation and interaction lengths, normal vector
@@ -173,6 +186,9 @@ private:
void GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
float &absorberCorrection) const;
+ void GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
+ float &absorberCorrection) const;
/**
* @brief Get number of active layers from position of a calo hit to the edge of the detector
@@ -235,6 +251,8 @@ private:
std::string m_encoder_str_LCal ;
std::string m_encoder_str_LHCal;
gear::GearMgr* _GEAR;
+ IGeomSvc* m_geosvc;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
};
//------------------------------------------------------------------------------------------------------------------------------------------
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/GeometryCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/GeometryCreator.h
index bb43d8a0f86a5211485a1e4cb301e3886da941a9..de894a78c27d288fc54028afa002d2beedd9a1cb 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/GeometryCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/GeometryCreator.h
@@ -11,6 +11,11 @@
#include "Api/PandoraApi.h"
#include "GaudiKernel/ISvcLocator.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/DetType.h"
+#include "DD4hep/DetectorSelector.h"
namespace gear { class CalorimeterParameters; class GearMgr; }
//------------------------------------------------------------------------------------------------------------------------------------------
@@ -53,6 +58,7 @@ public:
float m_hCalRingInnerPhiCoordinate; ///< HCal ring inner phi coordinate (missing from ILD gear files)
int m_hCalRingOuterSymmetryOrder; ///< HCal ring outer symmetry order (missing from ILD gear files)
float m_hCalRingOuterPhiCoordinate; ///< HCal ring outer phi coordinate (missing from ILD gear files)
+ bool m_use_dd4hep_geo; ///true to use dd4hep geo info
};
/**
@@ -102,6 +108,8 @@ private:
void SetDefaultSubDetectorParameters(const gear::CalorimeterParameters &inputParameters, const std::string &subDetectorName,
const pandora::SubDetectorType subDetectorType, PandoraApi::Geometry::SubDetector::Parameters ¶meters) const;
+ void SetDefaultSubDetectorParameters(const dd4hep::rec::LayeredCalorimeterData &inputParameters, const std::string &subDetectorName,
+ const pandora::SubDetectorType subDetectorType, PandoraApi::Geometry::SubDetector::Parameters ¶meters) const;
/**
* @brief Set positions of gaps in ILD detector and add information missing from GEAR parameters file
*
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
index 8188fb965dab87b4a2cdefcae3be7c24d7bccea2..c3c2599e7d518082db1d1981f388341bef069e26 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/PandoraPFAlg.h
@@ -273,6 +273,9 @@ protected:
std::vector<float> m_mc_charge;
int m_hasConversion;
+ Gaudi::Property<bool> m_use_dd4hep_geo{this, "use_dd4hep_geo", false,"choose if use geo info from dd4hep"};
+ Gaudi::Property<bool> m_use_dd4hep_decoder {this, "use_dd4hep_decoder", true,"if use decoder from dd4hep"};
+ Gaudi::Property<bool> m_use_preshower {this, "use_preshower", false,"if use preshower layer for calorimeter"};
Gaudi::Property< std::string > m_AnaOutput{ this, "AnaOutput", "/junofs/users/wxfang/MyGit/CEPCSW/Reconstruction/PFA/Pandora/GaudiPandora/Ana.root" };
//######################
std::map< std::string, std::string > m_collections;
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
index 134ff85bc5f08b7643155624e4e82910e7cbdce7..2b2c4011d2b24434816cd6bc94050dd540864a5e 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/TrackCreator.h
@@ -104,6 +104,7 @@ public:
float m_maxTpcInnerRDistance; ///< Track cut on distance from tpc inner r to id whether track can form pfo
float m_minTpcHitFractionOfExpected; ///< Minimum fraction of TPC hits compared to expected
int m_minFtdHitsForTpcHitFraction; ///< Minimum number of FTD hits to ignore TPC hit fraction
+ bool m_use_dd4hep_geo; ///
};
/**
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h b/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h
index 0bd75b5d94f4725581e44176a02fbb94a99088fa..07ae48091b2be4ea02d2bd79b547afdf9f3f0912 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/include/Utility.h
@@ -2,5 +2,17 @@
#define MYUTILITY 1
#include <sstream>
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/DetType.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/DetectorSelector.h"
std::string Convert (float number);
+dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, unsigned int excludeFlag) ;
+void line_a_b(float x1, float y1, float x2, float y2, float& a, float& b);
+float getPhi(float x, float y);
+int partition(float x, float y);
+int getLayer(float x, float y, float z, std::vector<float>& layers);
+int getLayer_v1(float x, float y, float z, std::vector<float>& layers);
+int getStave(float x, float y);
#endif
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
index e1c24b1e182b76ab503e007a0ec2d402424fe117..02009a5bdc83e9930e4f7b3642ddcdded49b300b 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/CaloHitCreator.cpp
@@ -49,14 +49,25 @@ CaloHitCreator::CaloHitCreator(const Settings &settings, const pandora::Pandora
throw "Failed to find GearSvc ...";
}
_GEAR = iSvc->getGearMgr();
+ if(m_settings.m_use_dd4hep_geo){
+ sc = svcloc->service("GeomSvc", m_geosvc, false);
+ if (!sc) throw "Failed to find GeomSvc.";
+ const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL),
+ ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) );
+ m_eCalBarrelOuterZ = eCalBarrelExtension->extent[3]/dd4hep::mm;
+ m_eCalBarrelInnerPhi0 = eCalBarrelExtension->inner_phi0/dd4hep::rad;
+ m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry;
+ }
+ else{
+ m_eCalBarrelOuterZ = (_GEAR->getEcalBarrelParameters().getExtent()[3]);
+ m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0());
+ m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder());
+ }
- m_eCalBarrelOuterZ = (_GEAR->getEcalBarrelParameters().getExtent()[3]);
m_hCalBarrelOuterZ = (_GEAR->getHcalBarrelParameters().getExtent()[3]);
m_muonBarrelOuterZ = (_GEAR->getYokeBarrelParameters().getExtent()[3]);
m_coilOuterR = (_GEAR->getGearParameters("CoilParameters").getDoubleVal("Coil_cryostat_outer_radius"));
- m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0());
- m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder());
m_hCalBarrelInnerPhi0 = (_GEAR->getHcalBarrelParameters().getPhi0());
m_hCalBarrelInnerSymmetry = (_GEAR->getHcalBarrelParameters().getSymmetryOrder());
m_muonBarrelInnerPhi0 = (_GEAR->getYokeBarrelParameters().getPhi0());
@@ -129,6 +140,12 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col
const std::string layerCodingString(m_encoder_str);
const std::string layerCoding(this->GetLayerCoding(layerCodingString));
const std::string staveCoding(this->GetStaveCoding(layerCodingString));
+ // get the DD4hep readout
+ const std::string name_readout = "EcalBarrelCollection";
+ if(m_settings.m_use_dd4hep_geo && m_settings.m_use_dd4hep_decoder ){
+ m_decoder = m_geosvc->getDecoder(name_readout);
+ if (!m_decoder) throw "Failed to get the decoder. ";
+ }
for (int i = 0; i < nElements; ++i)
{
@@ -173,7 +190,19 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col
PandoraApi::CaloHit::Parameters caloHitParameters;
caloHitParameters.m_hitType = pandora::ECAL;
caloHitParameters.m_isDigital = false;
- caloHitParameters.m_layer = cellIdDecoder(pCaloHit)[layerCoding.c_str()] + 1;
+ caloHitParameters.m_layer = m_settings.m_use_dd4hep_decoder == false ? cellIdDecoder(pCaloHit)[layerCoding.c_str()] + 1 : m_decoder->get(pCaloHit->getCellID(), "layer");// from 0 to 29, 0 is preshower layer
+ int Stave = 0 ;
+ if (m_settings.m_use_dd4hep_decoder == false){
+ Stave = cellIdDecoder(pCaloHit)[ staveCoding];
+ }
+ else{
+ Stave = m_decoder->get(pCaloHit->getCellID(), "stave");
+ Stave = Stave <=2 ? Stave+5 : Stave-3 ;//change to correct style
+ }
+ //std::cout<<"0Stave="<<Stave<<",0layer="<<caloHitParameters.m_layer.Get()<<std::endl;
+ if (Stave<0) throw "wrong Stave";
+ if (m_settings.m_use_preshower==false && caloHitParameters.m_layer.Get()<1) continue;//don't use preshower layer
+ //std::cout<<"Stave="<<Stave<<",layer="<<caloHitParameters.m_layer.Get()<<std::endl;
caloHitParameters.m_isInOuterSamplingLayer = false;
this->GetCommonCaloHitProperties(pCaloHit, caloHitParameters);
@@ -181,14 +210,27 @@ pandora::StatusCode CaloHitCreator::CreateECalCaloHits(const CollectionMaps& col
if (std::fabs(pCaloHit->getPosition()[2]) < m_eCalBarrelOuterZ)
{
- this->GetBarrelCaloHitProperties(pCaloHit, barrelLayerLayout, m_eCalBarrelInnerSymmetry, m_eCalBarrelInnerPhi0,
- cellIdDecoder(pCaloHit)[ staveCoding], caloHitParameters, absorberCorrection);
+ if(m_settings.m_use_dd4hep_geo){
+ const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer>& barrelLayers= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )->layers;
+ this->GetBarrelCaloHitProperties(pCaloHit, barrelLayers, m_eCalBarrelInnerSymmetry, m_eCalBarrelInnerPhi0,
+ Stave, caloHitParameters, absorberCorrection);
+ }
+ else{
+ this->GetBarrelCaloHitProperties(pCaloHit, barrelLayerLayout, m_eCalBarrelInnerSymmetry, m_eCalBarrelInnerPhi0,
+ Stave, caloHitParameters, absorberCorrection);
+ }
caloHitParameters.m_hadronicEnergy = eCalToHadGeVBarrel * pCaloHit->getEnergy();
}
else
{
- this->GetEndCapCaloHitProperties(pCaloHit, endcapLayerLayout, caloHitParameters, absorberCorrection);
+ if(m_settings.m_use_dd4hep_geo){
+ const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer>& endcapLayers= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )->layers;
+ this->GetEndCapCaloHitProperties(pCaloHit, endcapLayers, caloHitParameters, absorberCorrection);
+ }
+ else{
+ this->GetEndCapCaloHitProperties(pCaloHit, endcapLayerLayout, caloHitParameters, absorberCorrection);
+ }
caloHitParameters.m_hadronicEnergy = eCalToHadGeVEndCap * pCaloHit->getEnergy();
}
@@ -623,6 +665,57 @@ void CaloHitCreator::GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *c
pandora::CartesianVector(0, 0, -1);
}
+//------------------------------------------------------------------------------------------------------------------------------------------
+void CaloHitCreator::GetEndCapCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const
+{
+ caloHitParameters.m_hitRegion = pandora::ENDCAP;
+
+ const int physicalLayer(std::min(static_cast<int>(caloHitParameters.m_layer.Get()), static_cast<int>(layers.size()-1)));
+ caloHitParameters.m_cellSize0 = layers[physicalLayer].cellSize0/dd4hep::mm;
+ caloHitParameters.m_cellSize1 = layers[physicalLayer].cellSize1/dd4hep::mm;
+ //std::cout<<"Endcap m_cellSize0="<<caloHitParameters.m_cellSize0.Get()<<",m_cellSize1="<<caloHitParameters.m_cellSize1.Get()<<std::endl;
+ double thickness = (layers[physicalLayer].inner_thickness+layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ double nRadLengths = layers[physicalLayer].inner_nRadiationLengths;
+ double nIntLengths = layers[physicalLayer].inner_nInteractionLengths;
+ double layerAbsorberThickness = (layers[physicalLayer].inner_thickness-layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+
+ if(physicalLayer>0){
+ thickness += (layers[physicalLayer-1].outer_thickness -layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ nRadLengths += layers[physicalLayer-1].outer_nRadiationLengths;
+ nIntLengths += layers[physicalLayer-1].outer_nInteractionLengths;
+ layerAbsorberThickness += (layers[physicalLayer-1].outer_thickness -layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+
+ }
+ caloHitParameters.m_cellThickness = thickness;
+ caloHitParameters.m_nCellRadiationLengths = nRadLengths;
+ caloHitParameters.m_nCellInteractionLengths = nIntLengths;
+ if (caloHitParameters.m_nCellRadiationLengths.Get() < std::numeric_limits<float>::epsilon() || caloHitParameters.m_nCellInteractionLengths.Get() < std::numeric_limits<float>::epsilon())
+ {
+ std::cout<<"WARNING CaloHitCreator::GetEndCapCaloHitProperties Calo hit has 0 radiation length or interaction length: \
+ not creating a Pandora calo hit." << std::endl;
+ throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
+ }
+
+ absorberCorrection = 1.;
+ for (unsigned int i = 0, iMax = layers.size(); i < iMax; ++i)
+ {
+ float absorberThickness((layers[i].inner_thickness - layers[i].sensitive_thickness/2.0 )/dd4hep::mm);
+
+ if (i>0)
+ absorberThickness += (layers[i-1].outer_thickness - layers[i-1].sensitive_thickness/2.0)/dd4hep::mm;
+
+ if (absorberThickness < std::numeric_limits<float>::epsilon())
+ continue;
+
+ if (layerAbsorberThickness > std::numeric_limits<float>::epsilon())
+ absorberCorrection = absorberThickness / layerAbsorberThickness;
+
+ break;
+ }
+ caloHitParameters.m_cellNormalVector = (pCaloHit->getPosition()[2] > 0) ? pandora::CartesianVector(0, 0, 1) :
+ pandora::CartesianVector(0, 0, -1);
+}
//------------------------------------------------------------------------------------------------------------------------------------------
void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const gear::LayerLayout &layerLayout,
@@ -688,6 +781,81 @@ void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *c
}
}
+void CaloHitCreator::GetBarrelCaloHitProperties(const edm4hep::CalorimeterHit *const pCaloHit, const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer> &layers,
+ unsigned int barrelSymmetryOrder, float barrelPhi0, unsigned int staveNumber, PandoraApi::CaloHit::Parameters &caloHitParameters,
+ float &absorberCorrection) const
+{
+ caloHitParameters.m_hitRegion = pandora::BARREL;
+ const int physicalLayer(std::min(static_cast<int>(caloHitParameters.m_layer.Get()), static_cast<int>(layers.size()-1)));
+ caloHitParameters.m_cellSize0 = layers[physicalLayer].cellSize0/dd4hep::mm;
+ caloHitParameters.m_cellSize1 = layers[physicalLayer].cellSize1/dd4hep::mm;
+ std::cout<<"DD m_cellSize0="<<caloHitParameters.m_cellSize0.Get()<<",m_cellSize1="<<caloHitParameters.m_cellSize1.Get()<<std::endl;
+ double thickness = (layers[physicalLayer].inner_thickness+layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ double nRadLengths = layers[physicalLayer].inner_nRadiationLengths;
+ double nIntLengths = layers[physicalLayer].inner_nInteractionLengths;
+
+ double layerAbsorberThickness = (layers[physicalLayer].inner_thickness-layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ if(physicalLayer>0){
+ thickness += (layers[physicalLayer-1].outer_thickness -layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ nRadLengths += layers[physicalLayer-1].outer_nRadiationLengths;
+ nIntLengths += layers[physicalLayer-1].outer_nInteractionLengths;
+ layerAbsorberThickness += (layers[physicalLayer-1].outer_thickness -layers[physicalLayer].sensitive_thickness/2.0)/dd4hep::mm;
+ }
+
+ caloHitParameters.m_cellThickness = thickness;
+ caloHitParameters.m_nCellRadiationLengths = nRadLengths;
+ caloHitParameters.m_nCellInteractionLengths = nIntLengths;
+
+ if (caloHitParameters.m_nCellRadiationLengths.Get() < std::numeric_limits<float>::epsilon() || caloHitParameters.m_nCellInteractionLengths.Get() < std::numeric_limits<float>::epsilon())
+ {
+ std::cout<<"WARNIN CaloHitCreator::GetBarrelCaloHitProperties Calo hit has 0 radiation length or interaction length: \
+ not creating a Pandora calo hit." << std::endl;
+ throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER);
+ }
+ absorberCorrection = 1.;
+ float absorberThickness_0 = 0;
+ for (unsigned int i = 0, iMax = layers.size(); i < iMax; ++i)
+ {
+ float absorberThickness((layers[i].inner_thickness - layers[i].sensitive_thickness/2.0 )/dd4hep::mm);
+
+ if (i>0)
+ absorberThickness += (layers[i-1].outer_thickness - layers[i-1].sensitive_thickness/2.0)/dd4hep::mm;
+
+ if (absorberThickness < std::numeric_limits<float>::epsilon())
+ continue;
+
+ if (layerAbsorberThickness > std::numeric_limits<float>::epsilon())
+ absorberCorrection = absorberThickness / layerAbsorberThickness;
+ absorberThickness_0 = absorberThickness;
+ break;
+ }
+
+
+ //std::cout<<"DD m_cellSize0="<<caloHitParameters.m_cellSize0.Get()<<",m_cellSize1="<<caloHitParameters.m_cellSize1.Get()<<",physicalLayer="<<physicalLayer<<",layerAbsorberThickness="<<layerAbsorberThickness<<",absorberThickness_0="<<absorberThickness_0<<",m_cellThickness="<<caloHitParameters.m_cellThickness.Get()<<",m_nCellRadiationLengths="<<caloHitParameters.m_nCellRadiationLengths.Get()<<",m_nCellInteractionLengths="<<caloHitParameters.m_nCellInteractionLengths.Get()<<",absorberCorrection="<<absorberCorrection<<std::endl;
+
+
+ if (barrelSymmetryOrder > 2)
+ {
+ const float phi = barrelPhi0 + (2. * M_PI * static_cast<float>(staveNumber) / static_cast<float>(barrelSymmetryOrder));
+ caloHitParameters.m_cellNormalVector = pandora::CartesianVector(-std::sin(phi), std::cos(phi), 0);
+ }
+ else
+ {
+ const float pCaloHitPosition[3]={pCaloHit->getPosition()[0], pCaloHit->getPosition()[1], pCaloHit->getPosition()[2]};
+
+ if (pCaloHitPosition[1] != 0)
+ {
+ const float phi = barrelPhi0 + std::atan(pCaloHitPosition[0] / pCaloHitPosition[1]);
+ caloHitParameters.m_cellNormalVector = pandora::CartesianVector(std::sin(phi), std::cos(phi), 0);
+ }
+ else
+ {
+ caloHitParameters.m_cellNormalVector = (pCaloHitPosition[0] > 0) ? pandora::CartesianVector(1, 0, 0) :
+ pandora::CartesianVector(-1, 0, 0);
+ }
+ }
+}
+
//------------------------------------------------------------------------------------------------------------------------------------------
int CaloHitCreator::GetNLayersFromEdge(const edm4hep::CalorimeterHit *const pCaloHit) const
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp
index 4c88f7ba826fcf806b9a25f253cf51e06d1a92e2..4c3bcb3673d6c531060de032c9f9b943f12ede30 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/GeometryCreator.cpp
@@ -14,6 +14,7 @@
#include "gear/PadRowLayout2D.h"
#include "gear/LayerLayout.h"
#include "GeometryCreator.h"
+#include "Utility.h"
GeometryCreator::GeometryCreator(const Settings &settings, const pandora::Pandora *const pPandora) :
m_settings(settings),
@@ -73,9 +74,16 @@ void GeometryCreator::SetMandatorySubDetectorParameters(SubDetectorTypeMap &subD
std::cout << "Begin SetMandatorySubDetectorParameters:" << std::endl;
PandoraApi::Geometry::SubDetector::Parameters eCalBarrelParameters, eCalEndCapParameters, hCalBarrelParameters, hCalEndCapParameters,
muonBarrelParameters, muonEndCapParameters;
+ if(m_settings.m_use_dd4hep_geo){
+ std::cout << "Use dd4hep geo info for ECAL" << std::endl;
+ this->SetDefaultSubDetectorParameters(*const_cast<dd4hep::rec::LayeredCalorimeterData*>(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "ECalBarrel", pandora::ECAL_BARREL, eCalBarrelParameters);
+ this->SetDefaultSubDetectorParameters(*const_cast<dd4hep::rec::LayeredCalorimeterData*>(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "ECalEndCap", pandora::ECAL_ENDCAP, eCalEndCapParameters);
+ }
+ else{
+ this->SetDefaultSubDetectorParameters(_GEAR->getEcalBarrelParameters(), "ECalBarrel", pandora::ECAL_BARREL, eCalBarrelParameters);
+ this->SetDefaultSubDetectorParameters(_GEAR->getEcalEndcapParameters(), "ECalEndCap", pandora::ECAL_ENDCAP, eCalEndCapParameters);
+ }
- this->SetDefaultSubDetectorParameters(_GEAR->getEcalBarrelParameters(), "ECalBarrel", pandora::ECAL_BARREL, eCalBarrelParameters);
- this->SetDefaultSubDetectorParameters(_GEAR->getEcalEndcapParameters(), "ECalEndCap", pandora::ECAL_ENDCAP, eCalEndCapParameters);
this->SetDefaultSubDetectorParameters(_GEAR->getHcalBarrelParameters(), "HCalBarrel", pandora::HCAL_BARREL, hCalBarrelParameters);
this->SetDefaultSubDetectorParameters(_GEAR->getHcalEndcapParameters(), "HCalEndCap", pandora::HCAL_ENDCAP, hCalEndCapParameters);
this->SetDefaultSubDetectorParameters(_GEAR->getYokeBarrelParameters(), "MuonBarrel", pandora::MUON_BARREL, muonBarrelParameters);
@@ -209,6 +217,48 @@ void GeometryCreator::SetDefaultSubDetectorParameters(const gear::CalorimeterPar
}
}
+void GeometryCreator::SetDefaultSubDetectorParameters(const dd4hep::rec::LayeredCalorimeterData &inputParameters, const std::string &subDetectorName,
+ const pandora::SubDetectorType subDetectorType, PandoraApi::Geometry::SubDetector::Parameters ¶meters) const
+{
+ const std::vector<dd4hep::rec::LayeredCalorimeterStruct::Layer>& layers= inputParameters.layers;
+ parameters.m_subDetectorName = subDetectorName;
+ parameters.m_subDetectorType = subDetectorType;
+ parameters.m_innerRCoordinate = inputParameters.extent[0]/dd4hep::mm;
+ parameters.m_innerZCoordinate = inputParameters.extent[2]/dd4hep::mm;
+ parameters.m_innerPhiCoordinate = inputParameters.inner_phi0/dd4hep::rad;
+ parameters.m_innerSymmetryOrder = inputParameters.inner_symmetry;
+ parameters.m_outerRCoordinate = inputParameters.extent[1]/dd4hep::mm;
+ parameters.m_outerZCoordinate = inputParameters.extent[3]/dd4hep::mm;
+ parameters.m_outerPhiCoordinate = inputParameters.outer_phi0/dd4hep::rad;
+ parameters.m_outerSymmetryOrder = inputParameters.outer_symmetry;
+ parameters.m_isMirroredInZ = true;
+ parameters.m_nLayers = layers.size();
+
+ std::cout<<"DD m_subDetectorName="<<subDetectorName<<",m_subDetectorType="<<subDetectorType<<",m_innerRCoordinate="<<parameters.m_innerRCoordinate.Get()<<",m_innerZCoordinate="<<parameters.m_innerZCoordinate.Get()<<",m_innerPhiCoordinate="<<parameters.m_innerPhiCoordinate.Get()<<",m_innerSymmetryOrder="<<parameters.m_innerSymmetryOrder.Get()<<",m_outerRCoordinate="<<parameters.m_outerRCoordinate.Get()<<",m_outerZCoordinate="<<parameters.m_outerZCoordinate.Get()<<",m_outerPhiCoordinate="<<parameters.m_outerPhiCoordinate.Get()<<",m_outerSymmetryOrder="<<parameters.m_outerSymmetryOrder.Get()<<",m_nLayers="<<parameters.m_nLayers.Get()<<std::endl;
+ for (size_t i = 0; i< layers.size(); i++)
+ {
+ const dd4hep::rec::LayeredCalorimeterStruct::Layer & theLayer = layers.at(i);
+
+ PandoraApi::Geometry::LayerParameters layerParameters;
+
+ double totalNumberOfRadLengths = theLayer.inner_nRadiationLengths;
+ double totalNumberOfIntLengths = theLayer.inner_nInteractionLengths;
+
+ if(i>0){
+ //Add the numbers from previous layer's outer side
+ totalNumberOfRadLengths += layers.at(i-1).outer_nRadiationLengths;
+ totalNumberOfIntLengths += layers.at(i-1).outer_nInteractionLengths;
+
+ }
+
+ layerParameters.m_closestDistanceToIp = (theLayer.distance+theLayer.inner_thickness)/dd4hep::mm; //Distance to center of sensitive element
+ layerParameters.m_nRadiationLengths = totalNumberOfRadLengths;
+ layerParameters.m_nInteractionLengths = totalNumberOfIntLengths;
+
+ parameters.m_layerParametersVector.push_back(layerParameters);
+ std::cout<<"layer="<<i<<",m_closestDistanceToIp="<<layerParameters.m_closestDistanceToIp.Get()<<",m_nRadiationLengths="<<layerParameters.m_nRadiationLengths.Get()<<",m_nInteractionLengths="<<layerParameters.m_nInteractionLengths.Get()<<",outer_thickness="<<theLayer.outer_thickness<<",inner_thickness="<<theLayer.inner_thickness<<",sensitive_thickness="<<theLayer.sensitive_thickness<<std::endl;
+ }
+}
//------------------------------------------------------------------------------------------------------------------------------------------
pandora::StatusCode GeometryCreator::SetILDSpecificGeometry(SubDetectorTypeMap &subDetectorTypeMap, SubDetectorNameMap &subDetectorNameMap) const
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
index 14d4ee1a0e82b492706e30bba8cce509823bbbe0..6d458c6064b52536bdf0208c66167e1cdf11c00b 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/PandoraPFAlg.cpp
@@ -156,12 +156,16 @@ StatusCode PandoraPFAlg::initialize()
m_trackCreatorSettings.m_prongVertexCollections = m_ProngVertexCollections;
m_trackCreatorSettings.m_splitVertexCollections = m_SplitVertexCollections;
m_trackCreatorSettings.m_v0VertexCollections = m_V0VertexCollections;
+ m_trackCreatorSettings.m_use_dd4hep_geo = m_use_dd4hep_geo;
m_caloHitCreatorSettings.m_eCalCaloHitCollections = m_ECalCaloHitCollections;
m_caloHitCreatorSettings.m_hCalCaloHitCollections = m_HCalCaloHitCollections;
m_caloHitCreatorSettings.m_lCalCaloHitCollections = m_LCalCaloHitCollections;
m_caloHitCreatorSettings.m_lHCalCaloHitCollections = m_LHCalCaloHitCollections;
m_caloHitCreatorSettings.m_muonCaloHitCollections = m_MuonCaloHitCollections;
+ m_caloHitCreatorSettings.m_use_dd4hep_geo = m_use_dd4hep_geo;
+ m_caloHitCreatorSettings.m_use_dd4hep_decoder = m_use_dd4hep_decoder ;
+ m_caloHitCreatorSettings.m_use_preshower = m_use_preshower ;
m_mcParticleCreatorSettings.m_mcParticleCollections = m_MCParticleCollections;
m_mcParticleCreatorSettings.m_CaloHitRelationCollections = m_RelCaloHitCollections;
m_mcParticleCreatorSettings.m_TrackRelationCollections = m_RelTrackCollections;
@@ -252,6 +256,7 @@ StatusCode PandoraPFAlg::initialize()
m_geometryCreatorSettings.m_hCalRingInnerPhiCoordinate = m_HCalRingInnerPhiCoordinate;
m_geometryCreatorSettings.m_hCalRingOuterSymmetryOrder = m_HCalRingOuterSymmetryOrder;
m_geometryCreatorSettings.m_hCalRingOuterPhiCoordinate = m_HCalRingOuterPhiCoordinate;
+ m_geometryCreatorSettings.m_use_dd4hep_geo = m_use_dd4hep_geo;
// For Strip Splitting method and also for hybrid ECAL
m_caloHitCreatorSettings.m_stripSplittingOn = m_StripSplittingOn;
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
index c3b83199123447b624bbafbe38ec0db9d2561080..c2c0757fd3f9df1c0eb1684d73a07c41a2d3430c 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/TrackCreator.cpp
@@ -18,6 +18,13 @@
#include "gear/FTDParameters.h"
#include "gear/FTDLayerLayout.h"
+#include "DD4hep/Detector.h"
+#include "DD4hep/DD4hepUnits.h"
+#include "DD4hep/DetType.h"
+#include "DDRec/DetectorData.h"
+#include "DD4hep/DetectorSelector.h"
+#include "Utility.h"
+
#include "GaudiKernel/IService.h"
#include "GearSvc/IGearSvc.h"
#include "PandoraPFAlg.h"
@@ -48,11 +55,23 @@ TrackCreator::TrackCreator(const Settings &settings, const pandora::Pandora *con
m_tpcOuterR = (_GEAR->getTPCParameters().getPadLayout().getPlaneExtent()[1]);
m_tpcMaxRow = (_GEAR->getTPCParameters().getPadLayout().getNRows());
m_tpcZmax = (_GEAR->getTPCParameters().getMaxDriftLength());
- m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder());
- m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0());
- m_eCalBarrelInnerR = (_GEAR->getEcalBarrelParameters().getExtent()[0]);
- m_eCalEndCapInnerZ = (_GEAR->getEcalEndcapParameters().getExtent()[2]);
// fg: FTD description in GEAR has changed ...
+ if(m_settings.m_use_dd4hep_geo){
+ const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL),
+ ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) );
+ const dd4hep::rec::LayeredCalorimeterData * eCalEndcapExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP),
+ ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) );
+ m_eCalBarrelInnerPhi0 = eCalBarrelExtension->inner_phi0/dd4hep::rad;
+ m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry;
+ m_eCalBarrelInnerR = eCalBarrelExtension->extent[0]/dd4hep::mm;
+ m_eCalEndCapInnerZ = eCalEndcapExtension->extent[2]/dd4hep::mm;
+ }
+ else{
+ m_eCalBarrelInnerSymmetry = (_GEAR->getEcalBarrelParameters().getSymmetryOrder());
+ m_eCalBarrelInnerPhi0 = (_GEAR->getEcalBarrelParameters().getPhi0());
+ m_eCalBarrelInnerR = (_GEAR->getEcalBarrelParameters().getExtent()[0]);
+ m_eCalEndCapInnerZ = (_GEAR->getEcalEndcapParameters().getExtent()[2]);
+ }
try
{
m_ftdInnerRadii = _GEAR->getGearParameters("FTD").getDoubleVals("FTDInnerRadius");
diff --git a/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp b/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp
index 819fa15d4771c80d65e6437d22ebd46120d311de..697681ff616c9f9cbaecaf95a38c52b9e4aeaa0c 100644
--- a/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp
+++ b/Reconstruction/PFA/Pandora/GaudiPandora/src/Utility.cpp
@@ -5,3 +5,153 @@ std::string Convert (float number){
buff<<number;
return buff.str();
}
+dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, unsigned int excludeFlag=0) {
+
+
+ dd4hep::rec::LayeredCalorimeterData * theExtension = 0;
+
+ dd4hep::Detector & mainDetector = dd4hep::Detector::getInstance();
+ const std::vector< dd4hep::DetElement>& theDetectors = dd4hep::DetectorSelector(mainDetector).detectors( includeFlag, excludeFlag );
+
+
+// std::cout << " getExtension : includeFlag: " << dd4hep::DetType( includeFlag ) << " excludeFlag: " << dd4hep::DetType( excludeFlag )
+// << " found : " << theDetectors.size() << " - first det: " << theDetectors.at(0).name() << std::endl ;
+
+ if( theDetectors.size() != 1 ){
+
+ std::stringstream es ;
+ es << " getExtension: selection is not unique (or empty) includeFlag: " << dd4hep::DetType( includeFlag ) << " excludeFlag: " << dd4hep::DetType( excludeFlag )
+ << " --- found detectors : " ;
+ for( unsigned i=0, N= theDetectors.size(); i<N ; ++i ){
+ es << theDetectors.at(i).name() << ", " ;
+ }
+ throw std::runtime_error( es.str() ) ;
+ }
+
+ theExtension = theDetectors.at(0).extension<dd4hep::rec::LayeredCalorimeterData>();
+
+ return theExtension;
+}
+
+
+int partition(float x, float y)
+{
+ float a1, b1, a2, b2 ;
+ if (x>1850 && x < 2020){
+ line_a_b(1850, 750 , 2020, 600 , a1, b1);
+ line_a_b(1850, -1000, 2020, -850, a2, b2);
+ if ( (a1*x + b1*y) < 1 && (a2*x+b2*y) < 1 ) return 1 ;
+ }
+ if (y < 1850 && x < 2020){
+ line_a_b(750 , 1850,2020, 600, a1, b1);
+ line_a_b(1000, 1850,2020, 850, a2, b2);
+ if( (a1*x + b1*y) > 1 && (a2*x+b2*y) < 1 ) return 2;
+ }
+ if (y < 2020 && y > 1850){
+ line_a_b( -750, 1850 , -600, 2020 , a1, b1);
+ line_a_b( 1000, 1850 , 850 , 2020 , a2, b2);
+ if( (a1*x + b1*y) < 1 && (a2*x+b2*y) < 1 ) return 3;
+ }
+ if (y < 2020 && x > -1850){
+ line_a_b( -1850, 750 , -600 , 2020 , a1, b1);
+ line_a_b( -1850, 1000 , -850 , 2020 , a2, b2);
+ if( (a1*x + b1*y) > 1 && (a2*x+b2*y) < 1 ) return 4;
+ }
+ if (x > -2020 && x < -1850){
+ line_a_b( -1850, -750 , -2020, -600 , a1, b1);
+ line_a_b( -1850, 1000 , -2020, 850 , a2, b2);
+ if( (a1*x + b1*y) < 1 && (a2*x+b2*y) < 1 ) return 5;
+ }
+ if (x > -2020 && y > -1850){
+ line_a_b( -750 , -1850 , -2020, -600 , a1, b1);
+ line_a_b( -1000, -1850 , -2020, -850 , a2, b2);
+ if( (a1*x + b1*y) > 1 && (a2*x+b2*y) < 1) return 6;
+ }
+ if (y > -2020 && y < -1850){
+ line_a_b( 750 , -1850 , 600, -2020 , a1, b1);
+ line_a_b( -1000, -1850 , -850, -2020 , a2, b2);
+ if( (a1*x + b1*y) < 1 && (a2*x+b2*y) < 1) return 7;
+ }
+ if (y > -2020 && x < 1850){
+ line_a_b( 1850, -750 , 600, -2020 , a1, b1);
+ line_a_b( 1850, -1000 , 850, -2020 , a2, b2);
+ if( (a1*x + b1*y) > 1 && (a2*x+b2*y) < 1) return 8;
+ }
+ return -1;
+}
+
+void line_a_b(float x1, float y1, float x2, float y2, float& a, float& b)
+{
+ a = (y2-y1)/(x1*y2-x2*y1);
+ b = (x1-x2)/(x1*y2-x2*y1);
+}
+
+float getPhi(float x, float y)
+{
+ if (x==0 && y>0) return 90;
+ else if(x==0 && y<0) return 270;
+ else if(x==0 && y==0) return 0;
+ float phi = atan(y/x)*180/acos(-1);
+ if (x<0) phi = phi + 180;
+ else if (x>0 && y<0) phi = phi + 360;
+ return phi;
+}
+
+int getStave(float x, float y){
+
+int part = partition(x, y);
+int Stave = part>=3 ? (part-2) : (part+6);// correct S is from 1 to 8, S-1 is from 0-7
+return Stave-1;
+}
+
+int getLayer(float x, float y, float z, std::vector<float>& layers){
+if(abs(z) < 2400){
+ float phi = getPhi(x, y);
+ int part = partition(x, y);
+ float rotated = (part-1)*45;
+ float r = sqrt(x*x + y*y);
+ const float PI = acos(-1);
+ float new_x = r*cos((phi-rotated)*PI/180);
+ float new_y = r*sin((phi-rotated)*PI/180);
+ if(new_x <=layers.at(0)) return 0;
+ if(new_x >=layers.at(layers.size()-1)) return 29;
+ for(unsigned i=0; i<(layers.size()-1); i++){
+ if(new_x >layers.at(i) && new_x <=layers.at(i+1) ) return i+1;
+ }
+}
+else{
+ if(abs(z) <=layers.at(0)) return 0;
+ if(abs(z) >=layers.at(layers.size()-1)) return 29;
+ for(unsigned i=0; i<(layers.size()-1); i++){
+ if(abs(z) >layers.at(i) && abs(z) <=layers.at(i+1) ) return i+1;
+ }
+
+}
+return -1;
+}
+
+int getLayer_v1(float x, float y, float z, std::vector<float>& layers){
+if(abs(z) < 2400){
+ float phi = getPhi(x, y);
+ int part = partition(x, y);
+ float rotated = (part-1)*45;
+ float r = sqrt(x*x + y*y);
+ const float PI = acos(-1);
+ float new_x = r*cos((phi-rotated)*PI/180);
+ float new_y = r*sin((phi-rotated)*PI/180);
+ if(new_x <=layers.at(0)) return 1;
+ if(new_x >=layers.at(layers.size()-1)) return layers.size();
+ for(unsigned i=0; i<(layers.size()-1); i++){
+ if(new_x >layers.at(i) && new_x <=layers.at(i+1) ) return i+2;
+ }
+}
+else{
+ if(abs(z) <=layers.at(0)) return 1;
+ if(abs(z) >=layers.at(layers.size()-1)) return layers.size();
+ for(unsigned i=0; i<(layers.size()-1); i++){
+ if(abs(z) >layers.at(i) && abs(z) <=layers.at(i+1) ) return i+2;
+ }
+
+}
+return -1;
+}