diff --git a/include/DDCaloHitCreator.h b/include/DDCaloHitCreator.h index 7acd489..6283dd7 100644 --- a/include/DDCaloHitCreator.h +++ b/include/DDCaloHitCreator.h @@ -30,6 +30,9 @@ typedef std::vector CalorimeterHitVector; class DDCaloHitCreator { public: + // give access to the private memebers from DDCaloHitCreatorALLEGRO + friend class DDCaloHitCreatorALLEGRO; + typedef std::vector StringVector; typedef std::vector FloatVector; @@ -111,6 +114,10 @@ class DDCaloHitCreator float m_hCalBarrelOuterPhi0; ///< HCal barrel outer phi0 coordinate unsigned int m_hCalBarrelOuterSymmetry; ///< HCal barrel outer symmetry order + bool m_useSystemId; ///< flag whether to use systemId or not to identify origin of the CaloHit + int m_ecalBarrelSystemId; ///< systemId of ECal Barrel + int m_hcalBarrelSystemId; ///< systemId of HCal Barrel + public: FloatVector m_eCalBarrelNormalVector; FloatVector m_hCalBarrelNormalVector; diff --git a/include/DDCaloHitCreatorALLEGRO.h b/include/DDCaloHitCreatorALLEGRO.h new file mode 100644 index 0000000..1a3f378 --- /dev/null +++ b/include/DDCaloHitCreatorALLEGRO.h @@ -0,0 +1,111 @@ +/** + * @file DDMarlinPandora/include/DDCaloHitCreatorALLEGRO.h + * + * @brief Header file for the calo hit creator class. + * + * $Log: $ + */ + +#ifndef DDCALO_HIT_CREATORALLEGRO_H +#define DDCALO_HIT_CREATORALLEGRO_H 1 + +#include "EVENT/CalorimeterHit.h" +#include "EVENT/LCEvent.h" + +#include + +#include "Api/PandoraApi.h" + +#include +#include +#include + +#include "DDCaloHitCreator.h" + + +/** + * @brief DDCaloHitCreator class + */ +class DDCaloHitCreatorALLEGRO : public DDCaloHitCreator +{ +public: + typedef std::vector StringVector; + typedef std::vector FloatVector; + + /** + * @brief Constructor + * + * @param settings the creator settings + * @param pPandora address of the relevant pandora instance + */ + DDCaloHitCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora); + + /** + * @brief Destructor + */ + ~DDCaloHitCreatorALLEGRO(); + + /** + * @brief Create calo hits + * + * @param pLCEvent the lcio event + */ + pandora::StatusCode CreateCaloHits(const EVENT::LCEvent *const pLCEvent); + +private: + /** + * @brief Get common calo hit properties: position, parent address, input energy and time + * + * @param pCaloHit the lcio calorimeter hit + * @param caloHitParameters the calo hit parameters to populate + */ + void GetCommonCaloHitProperties(const EVENT::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const; + + /** + * @brief Get end cap specific calo hit properties: cell size, absorber radiation and interaction lengths, normal vector + * + * @param pCaloHit the lcio calorimeter hit + * @param layers the vector of layers from DDRec extensions + * @param caloHitParameters the calo hit parameters to populate + * @param absorberCorrection to receive the absorber thickness correction for the mip equivalent energy + */ + void GetEndCapCaloHitProperties(const EVENT::CalorimeterHit *const pCaloHit, const std::vector &layers, + PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const; + + /** + * @brief Get barrel specific calo hit properties: cell size, absorber radiation and interaction lengths, normal vector + * + * @param pCaloHit the lcio calorimeter hit + * @param layers the vector of layers from DDRec extensions + * @param barrelSymmetryOrder the barrel order of symmetry + * @param caloHitParameters the calo hit parameters to populate + * @param normalVector is the normalVector to the sensitive layers in local coordinates + * @param absorberCorrection to receive the absorber thickness correction for the mip equivalent energy + */ + void GetBarrelCaloHitProperties( const EVENT::CalorimeterHit *const pCaloHit, + const std::vector &layers, + unsigned int barrelSymmetryOrder, + PandoraApi::CaloHit::Parameters &caloHitParameters, + FloatVector const& normalVector, + float &absorberCorrection ) const; + + /** + * @brief Get number of active layers from position of a calo hit to the edge of the detector + * + * @param pCaloHit the lcio calorimeter hit + */ + int GetNLayersFromEdge(const EVENT::CalorimeterHit *const pCaloHit) const; + + /** + * @brief Get the maximum radius of a calo hit in a polygonal detector structure + * + * @param pCaloHit the lcio calorimeter hit + * @param symmetryOrder the symmetry order + * @param phi0 the angular orientation + * + * @return the maximum radius + */ + float GetMaximumRadius(const EVENT::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const; + +}; +#endif // #ifndef CALO_HIT_CREATOR_H diff --git a/include/DDGeometryCreator.h b/include/DDGeometryCreator.h index 733a611..9effbd4 100644 --- a/include/DDGeometryCreator.h +++ b/include/DDGeometryCreator.h @@ -21,6 +21,9 @@ class DDGeometryCreator { public: + // give access to the private memebers from DDGeometryCreatorALLEGRO + friend class DDGeometryCreatorALLEGRO; + /** * @brief Settings class */ diff --git a/include/DDGeometryCreatorALLEGRO.h b/include/DDGeometryCreatorALLEGRO.h new file mode 100644 index 0000000..e47035f --- /dev/null +++ b/include/DDGeometryCreatorALLEGRO.h @@ -0,0 +1,53 @@ +/** + * @file DDMarlinPandora/include/DDGeometryCreatorALLEGRO.h + * + * @brief Header file for the geometry creator class. + * + * $Log: $ + */ + +#ifndef DDGEOMETRYALLEGRO_CREATOR_H +#define DDGEOMETRYALLEGRO_CREATOR_H 1 + +#include "Api/PandoraApi.h" + +#include "DDRec/DetectorData.h" +#include "DDGeometryCreator.h" + +//------------------------------------------------------------------------------------------------------------------------------------------ + +/** + * @brief DDGeometryCreator class + */ +class DDGeometryCreatorALLEGRO : public DDGeometryCreator +{ +public: + /** + * @brief Constructor + * + * @param settings the creator settings + * @param pPandora address of the relevant pandora instance + */ + DDGeometryCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora); + + /** + * @brief Destructor + */ + ~DDGeometryCreatorALLEGRO(); + + /** + * @brief Create geometry + */ + pandora::StatusCode CreateGeometry() const; + +private: + /** + * @brief Set mandatory sub detector parameters + * + * @param subDetectorTypeMap the sub detector type map + */ + void SetMandatorySubDetectorParameters(SubDetectorTypeMap &subDetectorTypeMap) const; + +}; + +#endif // #ifndef GEOMETRY_CREATOR_H diff --git a/include/DDPandoraPFANewProcessor.h b/include/DDPandoraPFANewProcessor.h index 4f2d73d..ce9687c 100644 --- a/include/DDPandoraPFANewProcessor.h +++ b/include/DDPandoraPFANewProcessor.h @@ -11,6 +11,8 @@ #include "DDCaloHitCreator.h" #include "DDGeometryCreator.h" +#include "DDCaloHitCreatorALLEGRO.h" +#include "DDGeometryCreatorALLEGRO.h" #include "DDMCParticleCreator.h" #include "DDPfoCreator.h" #include "DDTrackCreatorBase.h" @@ -72,7 +74,8 @@ class DDPandoraPFANewProcessor : public marlin::Processor //Detector names not needed anymore, accessed by det type flags std::string m_trackCreatorName = ""; ///< The name of the DDTrackCreator implementation to use - + // Detector name (used by ALLEGRO) + std::string m_detectorName = ""; }; /** diff --git a/include/DDTrackCreatorALLEGRO.h b/include/DDTrackCreatorALLEGRO.h new file mode 100644 index 0000000..d0716c2 --- /dev/null +++ b/include/DDTrackCreatorALLEGRO.h @@ -0,0 +1,110 @@ +/** + * @file DDMarlinPandora/include/DDTrackCreatorALLEGRO.h + * + * @brief Header file for the ILD implementation of the track creator class. + * + * $Log: $ + */ + +#ifndef DDTRACK_CREATOR_ALLEGRO_H +#define DDTRACK_CREATOR_ALLEGRO_H 1 + +#include "DDTrackCreatorBase.h" + +//------------------------------------------------------------------------------------------------------------------------------------------ + +/** + * @brief DDTrackCreatorALLEGRO class + */ +class DDTrackCreatorALLEGRO : public DDTrackCreatorBase +{ +public: + + /** + * @brief Constructor + * + * @param settings the creator settings + * @param pPandora address of the relevant pandora instance + */ + DDTrackCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora); + + /** + * @brief Destructor + */ + virtual ~DDTrackCreatorALLEGRO(); + + /** + * @brief Create tracks implementation, insert user code here. Detector specific + * + * @param pLCEvent the lcio event + */ + pandora::StatusCode CreateTracks(EVENT::LCEvent *pLCEvent); + + + +protected: + + + //Detector-specific configuration variables + float m_trackerInnerR; ///< Inner radius of the barrel tracker + float m_trackerOuterR; ///< Outer radius of the barrel tracker + float m_trackerZmax; ///< max extent of the tracker along z + float m_cosTracker; + + DoubleVector m_endcapDiskInnerRadii; ///< List of endcapDisk inner radii + DoubleVector m_endcapDiskOuterRadii; ///< List of endcapDisk outer radii + DoubleVector m_endcapDiskZPositions; ///< List of endcapDisk z positions + unsigned int m_nEndcapDiskLayers; ///< Number of endcapDisk layers + unsigned int m_barrelTrackerLayers; ///< Number of barrel tracker layers + + float m_tanLambdaEndcapDisk; ///< Tan lambda for first endcapDisk layer + + /** + * @brief Whether track passes the quality cuts required in order to be used to form a pfo. Detector specific + * + * @param pTrack the lcio track + * @param trackParameters the track parameters + * + * @return boolean + */ + + virtual bool PassesQualityCuts(const EVENT::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const; + + /** + * @brief Decide whether track reaches the ecal surface. Detector specific + * + * @param pTrack the lcio track + * @param trackParameters the track parameters + */ + virtual void TrackReachesECAL(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; + + /** + * @brief Determine whether a track can be used to form a pfo under the following conditions: + * 1) if the track proves to be associated with a cluster, OR + * 2) if the track proves to have no cluster associations + * Detector specific + * + * @param pTrack the lcio track + * @param trackParameters the track parameters + */ + virtual void DefineTrackPfoUsage(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; + + /** + * @brief Copy track states stored in lcio tracks to pandora track parameters + * + * @param pTrack the lcio track + * @param trackParameters the track parameters + */ +// void GetTrackStates(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const; + + /** + * @brief Obtain track time when it reaches ECAL + * + * @param pTrack the lcio track + */ +// float CalculateTrackTimeAtCalorimeter(const EVENT::Track *const pTrack) const; +}; + + + +#endif // #ifndef DDTRACK_CREATOR_ALLEGRO_H diff --git a/src/DDCaloHitCreator.cc b/src/DDCaloHitCreator.cc index d29f180..51cf262 100644 --- a/src/DDCaloHitCreator.cc +++ b/src/DDCaloHitCreator.cc @@ -154,7 +154,8 @@ pandora::StatusCode DDCaloHitCreator::CreateECalCaloHits(const EVENT::LCEvent *c float absorberCorrection(1.); //FIXME: Why is this used to get barrel or endcap??? use SystemID if available - if (std::fabs(pCaloHit->getPosition()[2]) < m_settings.m_eCalBarrelOuterZ) + if ( (!m_settings.m_useSystemId && std::fabs(pCaloHit->getPosition()[2]) < m_settings.m_eCalBarrelOuterZ) || + (m_settings.m_useSystemId && cellIdDecoder(pCaloHit)["system"] == m_settings.m_ecalBarrelSystemId) ) { streamlog_out( DEBUG6 ) << "IDS " << *iter @@ -265,7 +266,8 @@ pandora::StatusCode DDCaloHitCreator::CreateHCalCaloHits(const EVENT::LCEvent *c float absorberCorrection(1.); - if (std::fabs(pCaloHit->getPosition()[2]) < m_settings.m_hCalBarrelOuterZ) + if ( (!m_settings.m_useSystemId && std::fabs(pCaloHit->getPosition()[2]) < m_settings.m_hCalBarrelOuterZ) || + (m_settings.m_useSystemId && cellIdDecoder(pCaloHit)["system"] == m_settings.m_hcalBarrelSystemId) ) { this->GetBarrelCaloHitProperties(pCaloHit, barrelLayers, m_settings.m_hCalBarrelInnerSymmetry, caloHitParameters, m_settings.m_hCalBarrelNormalVector, absorberCorrection); @@ -858,6 +860,9 @@ DDCaloHitCreator::Settings::Settings() m_hCalBarrelOuterR(0.f), m_hCalBarrelOuterPhi0(0.f), m_hCalBarrelOuterSymmetry(0.f), + m_useSystemId(false), + m_ecalBarrelSystemId(-1), + m_hcalBarrelSystemId(-1), m_eCalBarrelNormalVector({0.0, 0.0, 1.0}), m_hCalBarrelNormalVector({0.0, 0.0, 1.0}), m_muonBarrelNormalVector({0.0, 0.0, 1.0}) diff --git a/src/DDCaloHitCreatorALLEGRO.cc b/src/DDCaloHitCreatorALLEGRO.cc new file mode 100644 index 0000000..d59bb07 --- /dev/null +++ b/src/DDCaloHitCreatorALLEGRO.cc @@ -0,0 +1,295 @@ +/** + * @file DDMarlinPandora/src/DDCaloHitCreator.cc + * + * @brief Implementation of the calo hit creator class. + * + * $Log: $ + */ + +#include "DDCaloHitCreatorALLEGRO.h" + +#include "marlin/Global.h" +#include "marlin/Processor.h" + +#include "UTIL/CellIDDecoder.h" + + +#include +#include +#include +#include + +#include +#include +#include + +//forward declarations. See in DDPandoraPFANewProcessor.cc + +// dd4hep::rec::LayeredCalorimeterData * getExtension(std::string detectorName); +dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, unsigned int excludeFlag=0); + + +DDCaloHitCreatorALLEGRO::DDCaloHitCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora) : + DDCaloHitCreator(settings, pPandora) +{ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +DDCaloHitCreatorALLEGRO::~DDCaloHitCreatorALLEGRO() +{ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +pandora::StatusCode DDCaloHitCreatorALLEGRO::CreateCaloHits(const EVENT::LCEvent *const pLCEvent) +{ + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateECalCaloHits(pLCEvent)); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateHCalCaloHits(pLCEvent)); + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, this->CreateMuonCaloHits(pLCEvent)); + + return pandora::STATUS_CODE_SUCCESS; +} + + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDCaloHitCreatorALLEGRO::GetCommonCaloHitProperties(const EVENT::CalorimeterHit *const pCaloHit, PandoraApi::CaloHit::Parameters &caloHitParameters) const +{ + const float *pCaloHitPosition(pCaloHit->getPosition()); + const pandora::CartesianVector positionVector(pCaloHitPosition[0], pCaloHitPosition[1], pCaloHitPosition[2]); + + // FIXME! AD: for ECAL the cell gemoetry should be pandora::POINTING with cellSize0 = DeltaEta and cellSize1 = DeltaPhi + caloHitParameters.m_cellGeometry = pandora::RECTANGULAR; + caloHitParameters.m_positionVector = positionVector; + caloHitParameters.m_expectedDirection = positionVector.GetUnitVector(); + caloHitParameters.m_pParentAddress = (void*)pCaloHit; + caloHitParameters.m_inputEnergy = pCaloHit->getEnergy(); + caloHitParameters.m_time = pCaloHit->getTime(); +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDCaloHitCreatorALLEGRO::GetEndCapCaloHitProperties(const EVENT::CalorimeterHit *const pCaloHit, const std::vector &layers, + PandoraApi::CaloHit::Parameters &caloHitParameters, float &absorberCorrection) const +{ + caloHitParameters.m_hitRegion = pandora::ENDCAP; + + //FIXME! WHAT DO WE DO HERE? + const int physicalLayer(std::min(static_cast(caloHitParameters.m_layer.Get()), static_cast(layers.size()-1))); + caloHitParameters.m_cellSize0 = layers[physicalLayer].cellSize0/dd4hep::mm; + caloHitParameters.m_cellSize1 = layers[physicalLayer].cellSize1/dd4hep::mm; + + 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::epsilon() || caloHitParameters.m_nCellInteractionLengths.Get() < std::numeric_limits::epsilon()) + { + streamlog_out(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); + } + + //FIXME! do we need this? + 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::epsilon()) + continue; + + if (layerAbsorberThickness > std::numeric_limits::epsilon()) + absorberCorrection = absorberThickness / layerAbsorberThickness; + + break; + } + + caloHitParameters.m_cellNormalVector = (pCaloHit->getPosition()[2] > 0) ? pandora::CartesianVector(0, 0, 1) : + pandora::CartesianVector(0, 0, -1); +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDCaloHitCreatorALLEGRO::GetBarrelCaloHitProperties( const EVENT::CalorimeterHit *const pCaloHit, + const std::vector &layers, + unsigned int barrelSymmetryOrder, + PandoraApi::CaloHit::Parameters &caloHitParameters, + FloatVector const& normalVector, + float &absorberCorrection ) const +{ + caloHitParameters.m_hitRegion = pandora::BARREL; + + //FIXME! WHAT DO WE DO HERE? + const int physicalLayer(std::min(static_cast(caloHitParameters.m_layer.Get()), static_cast(layers.size()-1))); + caloHitParameters.m_cellSize0 = layers[physicalLayer].cellSize0/dd4hep::mm; + caloHitParameters.m_cellSize1 = layers[physicalLayer].cellSize1/dd4hep::mm; + + 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::epsilon() || caloHitParameters.m_nCellInteractionLengths.Get() < std::numeric_limits::epsilon()) + { + streamlog_out(WARNING) << "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); + } + + //FIXME! do we need this? + 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::epsilon()) + continue; + + if (layerAbsorberThickness > std::numeric_limits::epsilon()) + absorberCorrection = absorberThickness / layerAbsorberThickness; + + break; + } + + if (barrelSymmetryOrder > 2) + { + + if ( pCaloHit->getCellID0() != 0 ) { + + auto staveDetElement = m_volumeManager.lookupDetElement( pCaloHit->getCellID0() ); + dd4hep::Position local1(0.0, 0.0, 0.0); + dd4hep::Position local2(normalVector[0],normalVector[1],normalVector[2]); + dd4hep::Position global1(0.0, 0.0, 0.0); + dd4hep::Position global2(0.0, 0.0, 0.0); + staveDetElement.nominal().localToWorld( local1, global1 ); + staveDetElement.nominal().localToWorld( local2, global2 ); + dd4hep::Position normal( global2-global1 ); + + streamlog_out(DEBUG6) << " detelement: " << staveDetElement.name() + << " parent: " << staveDetElement.parent().name() + << " grandparent: " << staveDetElement.parent().parent().name() + << " cellID: " << pCaloHit->getCellID0() + << " PhiLoc:" << atan2( global1.y(), global1.x() )*180/M_PI + << " PhiNor:" << atan2( normal.y(), normal.x() )*180/M_PI + << " normal vector " + << std::setw(15) << normal.x() + << std::setw(15) << normal.y() + << std::setw(15) << normal.z() + << std::endl; + + caloHitParameters.m_cellNormalVector = pandora::CartesianVector( normal.x(), normal.y(), normal.z() ); + } else { + const double phi = atan2( pCaloHit->getPosition()[1], pCaloHit->getPosition()[0] ); + + streamlog_out( WARNING ) << "This hit does not have any cellIDs set, will use phi-direction for normal vector " + << " phi:" << std::setw(15) << phi*180/M_PI + << std::endl; + + caloHitParameters.m_cellNormalVector = pandora::CartesianVector( std::cos(phi), std::sin(phi) , 0.0 ); + } + + } + else + { + const float *pCaloHitPosition( pCaloHit->getPosition() ); + const float phi = std::atan2( pCaloHitPosition[1], pCaloHitPosition[0] ); + caloHitParameters.m_cellNormalVector = pandora::CartesianVector(std::cos(phi), std::sin(phi), 0); + } +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +int DDCaloHitCreatorALLEGRO::GetNLayersFromEdge(const EVENT::CalorimeterHit *const pCaloHit) const +{ + // Calo hit coordinate calculations + const float barrelMaximumRadius(this->GetMaximumRadius(pCaloHit, m_settings.m_hCalBarrelOuterSymmetry, m_settings.m_hCalBarrelOuterPhi0)); + const float endCapMaximumRadius(this->GetMaximumRadius(pCaloHit, m_settings.m_hCalEndCapInnerSymmetryOrder, m_settings.m_hCalEndCapInnerPhiCoordinate)); + const float caloHitAbsZ(std::fabs(pCaloHit->getPosition()[2])); + + // Distance from radial outer + float radialDistanceToEdge(std::numeric_limits::max()); + + if (caloHitAbsZ < m_settings.m_eCalBarrelOuterZ) + { + radialDistanceToEdge = (m_settings.m_hCalBarrelOuterR - barrelMaximumRadius) / m_hCalBarrelLayerThickness; + } + else + { + radialDistanceToEdge = (m_settings.m_hCalEndCapOuterR - endCapMaximumRadius) / m_hCalEndCapLayerThickness; + } + + // Distance from rear of endcap outer + float rearDistanceToEdge(std::numeric_limits::max()); + + if (caloHitAbsZ >= m_settings.m_eCalBarrelOuterZ) + { + rearDistanceToEdge = (m_settings.m_hCalEndCapOuterZ - caloHitAbsZ) / m_hCalEndCapLayerThickness; + } + else + { + const float rearDistance((m_settings.m_eCalBarrelOuterZ - caloHitAbsZ) / m_hCalBarrelLayerThickness); + + if (rearDistance < m_settings.m_layersFromEdgeMaxRearDistance) + { + const float overlapDistance((m_settings.m_hCalEndCapOuterR - endCapMaximumRadius) / m_hCalEndCapLayerThickness); + rearDistanceToEdge = std::max(rearDistance, overlapDistance); + } + } + + return static_cast(std::min(radialDistanceToEdge, rearDistanceToEdge)); +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +float DDCaloHitCreatorALLEGRO::GetMaximumRadius(const EVENT::CalorimeterHit *const pCaloHit, const unsigned int symmetryOrder, const float phi0) const +{ + const float *pCaloHitPosition(pCaloHit->getPosition()); + + if (symmetryOrder <= 2) + return std::sqrt((pCaloHitPosition[0] * pCaloHitPosition[0]) + (pCaloHitPosition[1] * pCaloHitPosition[1])); + + float maximumRadius(0.f); + const float twoPi(2.f * M_PI); + + for (unsigned int i = 0; i < symmetryOrder; ++i) + { + const float phi = phi0 + i * twoPi / static_cast(symmetryOrder); + float radius = pCaloHitPosition[0] * std::cos(phi) + pCaloHitPosition[1] * std::sin(phi); + + if (radius > maximumRadius) + maximumRadius = radius; + } + + return maximumRadius; +} diff --git a/src/DDGeometryCreatorALLEGRO.cc b/src/DDGeometryCreatorALLEGRO.cc new file mode 100644 index 0000000..f50e905 --- /dev/null +++ b/src/DDGeometryCreatorALLEGRO.cc @@ -0,0 +1,149 @@ +/** + * @file DDMarlinPandora/src/DDGeometryCreator.cc + * + * @brief Implementation of the geometry creator class. + * + * $Log: $ + */ + +#include "marlin/Global.h" +#include "marlin/Processor.h" + +#include "DDGeometryCreatorALLEGRO.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/DetectorData.h" +#include "DD4hep/DetType.h" +#include "DD4hep/DetectorSelector.h" + + +#include + +//Forward declarations. See DDPandoraPFANewProcessor.cc +// dd4hep::rec::LayeredCalorimeterData * getExtension(std::string detectorName); +dd4hep::rec::LayeredCalorimeterData * getExtension(unsigned int includeFlag, unsigned int excludeFlag=0); + +std::vector getTrackingRegionExtent(); + + +DDGeometryCreatorALLEGRO::DDGeometryCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora) : DDGeometryCreator(settings,pPandora) +{ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +DDGeometryCreatorALLEGRO::~DDGeometryCreatorALLEGRO() +{ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +pandora::StatusCode DDGeometryCreatorALLEGRO::CreateGeometry() const +{ + try + { + SubDetectorTypeMap subDetectorTypeMap; + this->SetMandatorySubDetectorParameters(subDetectorTypeMap); + + streamlog_out(DEBUG) << "Creating geometry for ALLEGRO detector"<< std::endl; + + for (SubDetectorTypeMap::const_iterator iter = subDetectorTypeMap.begin(), iterEnd = subDetectorTypeMap.end(); iter != iterEnd; ++iter) + PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::Geometry::SubDetector::Create(m_pPandora, iter->second)); + } + catch (std::exception &exception) + { + streamlog_out(ERROR) << "Failure in marlin pandora geometry creator, exception: " << exception.what() << std::endl; + throw exception; + } + + return pandora::STATUS_CODE_SUCCESS; +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDGeometryCreatorALLEGRO::SetMandatorySubDetectorParameters(SubDetectorTypeMap &subDetectorTypeMap) const +{ + PandoraApi::Geometry::SubDetector::Parameters eCalBarrelParameters, eCalEndCapParameters, hCalBarrelParameters, hCalEndCapParameters, + muonBarrelParameters, muonEndCapParameters; + + this->SetDefaultSubDetectorParameters(*const_cast(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(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "ECalEndCap", pandora::ECAL_ENDCAP, eCalEndCapParameters); + this->SetDefaultSubDetectorParameters(*const_cast(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::HADRONIC | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "HCalBarrel", pandora::HCAL_BARREL, hCalBarrelParameters); + this->SetDefaultSubDetectorParameters(*const_cast(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::HADRONIC | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "HCalEndCap", pandora::HCAL_ENDCAP, hCalEndCapParameters); + this->SetDefaultSubDetectorParameters(*const_cast(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::MUON | dd4hep::DetType::BARREL), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "MuonBarrel", pandora::MUON_BARREL, muonBarrelParameters); + this->SetDefaultSubDetectorParameters(*const_cast(getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::MUON | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) )), "MuonEndCap", pandora::MUON_ENDCAP, muonEndCapParameters); + + subDetectorTypeMap[pandora::ECAL_BARREL] = eCalBarrelParameters; + subDetectorTypeMap[pandora::ECAL_ENDCAP] = eCalEndCapParameters; + subDetectorTypeMap[pandora::HCAL_BARREL] = hCalBarrelParameters; + subDetectorTypeMap[pandora::HCAL_ENDCAP] = hCalEndCapParameters; + subDetectorTypeMap[pandora::MUON_BARREL] = muonBarrelParameters; + subDetectorTypeMap[pandora::MUON_ENDCAP] = muonEndCapParameters; + + // FIXME! AD: currently ignoring tracker parameters since we are using truth tracks + /* + PandoraApi::Geometry::SubDetector::Parameters trackerParameters; + + trackerParameters.m_subDetectorName = "Tracker"; + trackerParameters.m_subDetectorType = pandora::INNER_TRACKER; + trackerParameters.m_innerRCoordinate = getTrackingRegionExtent()[0]; + trackerParameters.m_innerZCoordinate = 0.f; + trackerParameters.m_innerPhiCoordinate = 0.f; + trackerParameters.m_innerSymmetryOrder = 0; + trackerParameters.m_outerRCoordinate = getTrackingRegionExtent()[1]; + trackerParameters.m_outerZCoordinate = getTrackingRegionExtent()[2]; + trackerParameters.m_outerPhiCoordinate = 0.f; + trackerParameters.m_outerSymmetryOrder = 0; + trackerParameters.m_isMirroredInZ = true; + trackerParameters.m_nLayers = 0; + subDetectorTypeMap[pandora::INNER_TRACKER] = trackerParameters; + */ + + ///FIXME:Implement a parameter for the Coil/Solenoid name + ///NOTE: Is this the way to go, or should we go with reco structure? + try{ + PandoraApi::Geometry::SubDetector::Parameters coilParameters; + + /* + const dd4hep::rec::LayeredCalorimeterData * coilExtension= getExtension( ( dd4hep::DetType::COIL ) ); + coilParameters.m_subDetectorName = "Coil"; + coilParameters.m_subDetectorType = pandora::COIL; + coilParameters.m_innerRCoordinate = coilExtension->extent[0]/ dd4hep::mm; + coilParameters.m_innerZCoordinate = 0.f; + coilParameters.m_innerPhiCoordinate = 0.f; + coilParameters.m_innerSymmetryOrder = 0; + coilParameters.m_outerRCoordinate = coilExtension->extent[1]/ dd4hep::mm; + coilParameters.m_outerZCoordinate = coilExtension->extent[3]/ dd4hep::mm; + coilParameters.m_outerPhiCoordinate = 0.f; + coilParameters.m_outerSymmetryOrder = 0; + coilParameters.m_isMirroredInZ = true; + coilParameters.m_nLayers = 0; + subDetectorTypeMap[pandora::COIL] = coilParameters; + */ + + // FIXME! AD: below an ugly way is used to define the geometry of the Solenoid since we do not have the LayeredCalorimeterData for Coil in ALLEGRO + //Get ECal Barrel extension + const dd4hep::rec::LayeredCalorimeterData * eCalBarrelExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::ELECTROMAGNETIC | dd4hep::DetType::BARREL), + ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); + //Get HCal Barrel extension + const dd4hep::rec::LayeredCalorimeterData * hCalBarrelExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::HADRONIC | dd4hep::DetType::BARREL), + ( dd4hep::DetType::AUXILIARY | dd4hep::DetType::FORWARD ) ); + coilParameters.m_subDetectorName = "Coil"; + coilParameters.m_subDetectorType = pandora::COIL; + coilParameters.m_innerRCoordinate = eCalBarrelExtension->extent[1]/ dd4hep::mm; // ECAL barrel outer R is assumed to be the inner R for the Solenoid + coilParameters.m_innerZCoordinate = 0.f; + coilParameters.m_innerPhiCoordinate = 0.f; + coilParameters.m_innerSymmetryOrder = 0; + coilParameters.m_outerRCoordinate = hCalBarrelExtension->extent[0]/ dd4hep::mm; // HCAL barrel inner R is assumed to be the outer R for the Solenoid + coilParameters.m_outerZCoordinate = eCalBarrelExtension->extent[3]/ dd4hep::mm; // Solenoid outer Z is assumed to be the same as ECAL barrel outer Z + coilParameters.m_outerPhiCoordinate = 0.f; + coilParameters.m_outerSymmetryOrder = 0; + coilParameters.m_isMirroredInZ = true; + coilParameters.m_nLayers = 0; + subDetectorTypeMap[pandora::COIL] = coilParameters; + + } catch ( std::exception & e ) { + streamlog_out(ERROR) << "Failed to access COIL parameters: "<FinaliseSteeringParameters(); m_pPandora = new pandora::Pandora(); - m_pGeometryCreator = new DDGeometryCreator(m_geometryCreatorSettings, m_pPandora); - m_pCaloHitCreator = new DDCaloHitCreator(m_caloHitCreatorSettings, m_pPandora); - - + if(m_settings.m_detectorName == "ALLEGRO") + { + m_settings.m_trackCreatorName = "DDTrackCreatorALLEGRO"; + m_pGeometryCreator = new DDGeometryCreatorALLEGRO(m_geometryCreatorSettings, m_pPandora); + m_pCaloHitCreator = new DDCaloHitCreatorALLEGRO(m_caloHitCreatorSettings, m_pPandora); + } + else + { + m_pGeometryCreator = new DDGeometryCreator(m_geometryCreatorSettings, m_pPandora); + m_pCaloHitCreator = new DDCaloHitCreator(m_caloHitCreatorSettings, m_pPandora); + } + ///FIXME: IMPLEMENT FACTORY if (m_settings.m_trackCreatorName == "DDTrackCreatorCLIC") m_pTrackCreator = new DDTrackCreatorCLIC(m_trackCreatorSettings, m_pPandora); else if (m_settings.m_trackCreatorName == "DDTrackCreatorILD") m_pTrackCreator = new DDTrackCreatorILD(m_trackCreatorSettings, m_pPandora); + else if (m_settings.m_trackCreatorName == "DDTrackCreatorALLEGRO") + m_pTrackCreator = new DDTrackCreatorALLEGRO(m_trackCreatorSettings, m_pPandora); else streamlog_out(ERROR) << "Unknown DDTrackCreator: "<RegisterUserComponents()); - PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pGeometryCreator->CreateGeometry()); + if(m_settings.m_detectorName == "ALLEGRO") + { + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, static_cast(m_pGeometryCreator)->CreateGeometry()); + } + else + { + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pGeometryCreator->CreateGeometry()); + } PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::ReadSettings(*m_pPandora, m_settings.m_pandoraSettingsXmlFile)); } catch (pandora::StatusCodeException &statusCodeException) @@ -212,10 +229,18 @@ void DDPandoraPFANewProcessor::processEvent(LCEvent *pLCEvent) (void) m_pandoraToLCEventMap.insert(PandoraToLCEventMap::value_type(m_pPandora, pLCEvent)); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pDDMCParticleCreator->CreateMCParticles(pLCEvent)); - PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTrackAssociations(pLCEvent)); - PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTracks(pLCEvent)); - PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pDDMCParticleCreator->CreateTrackToMCParticleRelationships(pLCEvent, m_pTrackCreator->GetTrackVector())); - PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pCaloHitCreator->CreateCaloHits(pLCEvent)); + if(m_settings.m_detectorName == "ALLEGRO") + { + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, static_cast(m_pTrackCreator)->CreateTracks(pLCEvent)); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, static_cast(m_pCaloHitCreator)->CreateCaloHits(pLCEvent)); + } + else + { + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTrackAssociations(pLCEvent)); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTracks(pLCEvent)); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pDDMCParticleCreator->CreateTrackToMCParticleRelationships(pLCEvent, m_pTrackCreator->GetTrackVector())); + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pCaloHitCreator->CreateCaloHits(pLCEvent)); + } PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pDDMCParticleCreator->CreateCaloHitToMCParticleRelationships(pLCEvent, m_pCaloHitCreator->GetCalorimeterHitVector())); PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::ProcessEvent(*m_pPandora)); @@ -816,6 +841,12 @@ void DDPandoraPFANewProcessor::ProcessSteeringFile() m_settings.m_trackCreatorName, std::string("DDTrackCreatorCLIC")); + ///EXTRA parameter that would initialize everything for ALLEGRO detector if m_detectorName=="ALLEGRO" + registerProcessorParameter("DetectorName", + "The name of the detector", + m_settings.m_detectorName, + std::string("")); + registerProcessorParameter("ECalBarrelNormalVector", "Normal vector for the ECal barrel sensitive layers in local coordinates", m_caloHitCreatorSettings.m_eCalBarrelNormalVector, @@ -900,19 +931,28 @@ void DDPandoraPFANewProcessor::FinaliseSteeringParameters() // //Get Muon Endcap extension by type, ignore plugs and rings // const dd4hep::rec::LayeredCalorimeterData * muonEndcapExtension= getExtension( ( dd4hep::DetType::CALORIMETER | dd4hep::DetType::MUON | dd4hep::DetType::ENDCAP), ( dd4hep::DetType::AUXILIARY ) ); - //Get COIL extension - const dd4hep::rec::LayeredCalorimeterData * coilExtension= getExtension( ( dd4hep::DetType::COIL ) ); - - + if(m_settings.m_detectorName != "ALLEGRO") + { + //Get COIL extension + const dd4hep::rec::LayeredCalorimeterData * coilExtension= getExtension( ( dd4hep::DetType::COIL ) ); + m_caloHitCreatorSettings.m_coilOuterR = coilExtension->extent[1]/dd4hep::mm; + } + else + { + m_caloHitCreatorSettings.m_useSystemId = true; + m_caloHitCreatorSettings.m_ecalBarrelSystemId = 4; + m_caloHitCreatorSettings.m_hcalBarrelSystemId = 8; + m_caloHitCreatorSettings.m_coilOuterR = 0; + } + m_trackCreatorSettings.m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry; m_trackCreatorSettings.m_eCalBarrelInnerPhi0 = eCalBarrelExtension->inner_phi0/dd4hep::rad; m_trackCreatorSettings.m_eCalBarrelInnerR = eCalBarrelExtension->extent[0]/dd4hep::mm; m_trackCreatorSettings.m_eCalEndCapInnerZ = eCalEndcapExtension->extent[2]/dd4hep::mm; - + m_caloHitCreatorSettings.m_eCalBarrelOuterZ = eCalBarrelExtension->extent[3]/dd4hep::mm; m_caloHitCreatorSettings.m_hCalBarrelOuterZ = hCalBarrelExtension->extent[3]/dd4hep::mm; m_caloHitCreatorSettings.m_muonBarrelOuterZ = muonBarrelExtension->extent[3]/dd4hep::mm; - m_caloHitCreatorSettings.m_coilOuterR = coilExtension->extent[1]/dd4hep::mm; m_caloHitCreatorSettings.m_eCalBarrelInnerPhi0 = eCalBarrelExtension->inner_phi0/dd4hep::rad; m_caloHitCreatorSettings.m_eCalBarrelInnerSymmetry = eCalBarrelExtension->inner_symmetry; m_caloHitCreatorSettings.m_hCalBarrelInnerPhi0 = hCalBarrelExtension->inner_phi0/dd4hep::rad; diff --git a/src/DDTrackCreatorALLEGRO.cc b/src/DDTrackCreatorALLEGRO.cc new file mode 100644 index 0000000..7d6c805 --- /dev/null +++ b/src/DDTrackCreatorALLEGRO.cc @@ -0,0 +1,425 @@ +/** + * @file DDMarlinPandora/src/DDTrackCreatorALLEGRO.cc + * + * @brief Implementation of the track creator class for a ALLEGRO all silicon tracker. + * + * $Log: $ + */ + +#include "DDTrackCreatorALLEGRO.h" + +#include "marlin/Global.h" +#include "marlin/Processor.h" + +#include "EVENT/LCCollection.h" +#include "EVENT/ReconstructedParticle.h" +#include "EVENT/Vertex.h" +#include "UTIL/ILDConf.h" +#include "UTIL/Operators.h" + +#include "Pandora/PdgTable.h" +#include "LCObjects/LCTrack.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/DD4hepUnits.h" +#include "DDRec/DetectorData.h" +#include "DD4hep/DetType.h" +#include "DD4hep/DetectorSelector.h" + +#include +#include +#include + +//forward declarations. See in DDPandoraPFANewProcessor.cc +std::vector getTrackingRegionExtent(); + +DDTrackCreatorALLEGRO::DDTrackCreatorALLEGRO(const Settings &settings, const pandora::Pandora *const pPandora) + : DDTrackCreatorBase(settings,pPandora), + m_trackerInnerR( 0.f ), + m_trackerOuterR( 0.f ), + m_trackerZmax( 0.f ), + m_cosTracker( 0.f ), + m_endcapDiskInnerRadii( DoubleVector() ), + m_endcapDiskOuterRadii( DoubleVector() ), + m_endcapDiskZPositions( DoubleVector() ), + m_nEndcapDiskLayers( 0 ), + m_barrelTrackerLayers( 0 ), + m_tanLambdaEndcapDisk( 0.f ) + +{ + // FIXME! AD: currently ignoring the tracker parameters since we are working with truth tracks... + /* + m_trackerInnerR = getTrackingRegionExtent()[0]; + m_trackerOuterR = getTrackingRegionExtent()[1]; + m_trackerZmax = getTrackingRegionExtent()[2]; + + ///FIXME: Probably need to be something relating to last disk inner radius + m_cosTracker = m_trackerZmax / std::sqrt(m_trackerZmax * m_trackerZmax + m_trackerInnerR * m_trackerInnerR); + + dd4hep::Detector & mainDetector = dd4hep::Detector::getInstance(); + + //Maybe we need to veto the vertex? That was done in the ILD case + const std::vector< dd4hep::DetElement>& barrelDets = dd4hep::DetectorSelector(mainDetector).detectors( ( dd4hep::DetType::TRACKER | dd4hep::DetType::BARREL )) ; + + m_barrelTrackerLayers = 0; + for (std::vector< dd4hep::DetElement>::const_iterator iter = barrelDets.begin(), iterEnd = barrelDets.end();iter != iterEnd; ++iter){ + try + { + dd4hep::rec::ZPlanarData * theExtension = 0; + + const dd4hep::DetElement& theDetector = *iter; + theExtension = theDetector.extension(); + + unsigned int N = theExtension->layers.size(); + m_barrelTrackerLayers=m_barrelTrackerLayers+N; + + streamlog_out( DEBUG2 ) << " Adding layers for barrel tracker from DD4hep for "<< theDetector.name()<< "- n layers: " << N<< " sum up to now: "<(*iter).name()<<" : " << exception.what() << std::endl; + } + } + + m_nEndcapDiskLayers=0; + m_endcapDiskInnerRadii.clear(); + m_endcapDiskOuterRadii.clear(); + + //Instead of gear, loop over a provided list of forward (read: endcap) tracking detectors. For ILD this would be FTD + const std::vector< dd4hep::DetElement>& endcapDets = dd4hep::DetectorSelector(mainDetector).detectors( ( dd4hep::DetType::TRACKER | dd4hep::DetType::ENDCAP )) ; + + for (std::vector< dd4hep::DetElement>::const_iterator iter = endcapDets.begin(), iterEnd = endcapDets.end();iter != iterEnd; ++iter){ + try + { + dd4hep::rec::ZDiskPetalsData * theExtension = 0; + + const dd4hep::DetElement& theDetector = *iter; + theExtension = theDetector.extension(); + + unsigned int N = theExtension->layers.size(); + streamlog_out( DEBUG2 ) << " Filling FTD-like parameters from DD4hep for "<< theDetector.name()<< "- n layers: " << N<< std::endl; + + for(unsigned int i = 0; i < N; ++i) + { + dd4hep::rec::ZDiskPetalsData::LayerLayout thisLayer = theExtension->layers[i]; + + // Create a disk to represent even number petals front side + //FIXME! VERIFY THAT TIS MAKES SENSE! + m_endcapDiskInnerRadii.push_back(thisLayer.distanceSensitive/dd4hep::mm); + m_endcapDiskOuterRadii.push_back(thisLayer.distanceSensitive/dd4hep::mm+thisLayer.lengthSensitive/dd4hep::mm); + + // Take the mean z position of the staggered petals + const double zpos(thisLayer.zPosition/dd4hep::mm); + m_endcapDiskZPositions.push_back(zpos); + streamlog_out( DEBUG2 ) << " layer " << i << " - mean z position = " << zpos << std::endl; + } + + m_nEndcapDiskLayers = m_endcapDiskZPositions.size() ; + } catch (std::runtime_error &exception){ + streamlog_out(WARNING) << "DDTrackCreatorALLEGRO exception during Forward Tracking Disk parameter construction for detector "<(*iter).name()<<" : " << exception.what() << std::endl; + } + } + + for (unsigned int iEndcapDiskLayer = 0; iEndcapDiskLayer < m_nEndcapDiskLayers; ++iEndcapDiskLayer) + { + if ((std::fabs(m_endcapDiskOuterRadii[iEndcapDiskLayer]) < std::numeric_limits::epsilon()) || + (std::fabs(m_endcapDiskInnerRadii[iEndcapDiskLayer]) < std::numeric_limits::epsilon())) + { + throw pandora::StatusCodeException(pandora::STATUS_CODE_INVALID_PARAMETER); + } + } + + m_tanLambdaEndcapDisk = m_endcapDiskZPositions[0] / m_endcapDiskOuterRadii[0]; + */ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +DDTrackCreatorALLEGRO::~DDTrackCreatorALLEGRO() +{ +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +pandora::StatusCode DDTrackCreatorALLEGRO::CreateTracks(EVENT::LCEvent *pLCEvent) +{ + for (StringVector::const_iterator iter = m_settings.m_trackCollections.begin(), iterEnd = m_settings.m_trackCollections.end(); + iter != iterEnd; ++iter) + { + try + { + const EVENT::LCCollection *pTrackCollection = pLCEvent->getCollection(*iter); + for (int i = 0, iMax = pTrackCollection->getNumberOfElements(); i < iMax; ++i) + { + EVENT::Track *pTrack = dynamic_cast(pTrackCollection->getElementAt(i)); + + if (NULL == pTrack) + throw EVENT::Exception("Collection type mismatch"); + + streamlog_out(DEBUG0)<<" Warning! Ignoring expected number of hits and other hit number cuts. Should eventually change!"<getD0(); + trackParameters.m_z0 = pTrack->getZ0(); + trackParameters.m_pParentAddress = pTrack; + + // By default, assume tracks are charged pions + const float signedCurvature(pTrack->getOmega()); + trackParameters.m_particleId = (signedCurvature > 0) ? pandora::PI_PLUS : pandora::PI_MINUS; + trackParameters.m_mass = pandora::PdgTable::GetParticleMass(pandora::PI_PLUS); + + // Use particle id information from V0 and Kink finders + TrackToPidMap::const_iterator trackPIDiter = m_trackToPidMap.find(pTrack); + + if(trackPIDiter != m_trackToPidMap.end()) + { + trackParameters.m_particleId = trackPIDiter->second; + trackParameters.m_mass = pandora::PdgTable::GetParticleMass(trackPIDiter->second); + } + + if (0.f != signedCurvature) + trackParameters.m_charge = static_cast(signedCurvature / std::fabs(signedCurvature)); + + try + { + this->GetTrackStates(pTrack, trackParameters); + this->TrackReachesECAL(pTrack, trackParameters); + // FIXME! AD: this is disabled for ALLEGRO since I am not sure what to do here with truth tracks + //this->GetTrackStatesAtCalo(pTrack, trackParameters); + this->DefineTrackPfoUsage(pTrack, trackParameters); + + PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::Track::Create(m_pandora, trackParameters, *m_lcTrackFactory)); + m_trackVector.push_back(pTrack); + } + catch (pandora::StatusCodeException &statusCodeException) + { + streamlog_out(ERROR) << "Failed to extract a track: " << statusCodeException.ToString() << std::endl; + streamlog_out( DEBUG3 ) << " failed track : " << *pTrack << std::endl ; + } + catch (EVENT::Exception &exception) + { + streamlog_out(WARNING) << "Failed to extract a vertex: " << exception.what() << std::endl; + } + } + streamlog_out( DEBUG5 ) << "After treating collection : " << *iter<<" with "<getNumberOfElements()<<" tracks, the track vector size is "<< m_trackVector.size()<< std::endl ; + + } + catch (EVENT::Exception &exception) + { + streamlog_out(WARNING) << "Failed to extract track collection: " << *iter << ", " << exception.what() << std::endl; + } + } + + return pandora::STATUS_CODE_SUCCESS; +} + +bool DDTrackCreatorALLEGRO::PassesQualityCuts(const EVENT::Track *const pTrack, const PandoraApi::Track::Parameters &trackParameters) const +{ + // First simple sanity checks + if (trackParameters.m_trackStateAtCalorimeter.Get().GetPosition().GetMagnitude() < m_settings.m_minTrackECalDistanceFromIp){ + streamlog_out(WARNING) << " Dropping track! Distance at ECAL: " << trackParameters.m_trackStateAtCalorimeter.Get().GetPosition().GetMagnitude()<getOmega() == 0.f) + { + streamlog_out(ERROR) << "Track has Omega = 0 " << std::endl; + return false; + } + + // Check momentum uncertainty is reasonable to use track + const pandora::CartesianVector &momentumAtDca(trackParameters.m_momentumAtDca.Get()); + const float sigmaPOverP(std::sqrt(pTrack->getCovMatrix()[5]) / std::fabs(pTrack->getOmega())); + + if (sigmaPOverP > m_settings.m_maxTrackSigmaPOverP) + { + streamlog_out(WARNING) << " Dropping track : " << momentumAtDca.GetMagnitude() << "+-" << sigmaPOverP * (momentumAtDca.GetMagnitude()) + << " chi2 = " << pTrack->getChi2() << " " << pTrack->getNdf() + << " from " << pTrack->getTrackerHits().size() << std::endl ; + + streamlog_out(DEBUG3) << " track : " << *pTrack + << std::endl; + return false; + } + + return true; +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDTrackCreatorALLEGRO::DefineTrackPfoUsage(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const +{ + bool canFormPfo(false); + bool canFormClusterlessPfo(false); + + if (trackParameters.m_reachesCalorimeter.Get() && !this->IsParent(pTrack)) + { + const float d0(std::fabs(pTrack->getD0())), z0(std::fabs(pTrack->getZ0())); + + EVENT::TrackerHitVec trackerHitvec(pTrack->getTrackerHits()); + float rInner(std::numeric_limits::max()), zMin(std::numeric_limits::max()); + + for (EVENT::TrackerHitVec::const_iterator iter = trackerHitvec.begin(), iterEnd = trackerHitvec.end(); iter != iterEnd; ++iter) + { + const double *pPosition((*iter)->getPosition()); + const float x(pPosition[0]), y(pPosition[1]), absoluteZ(std::fabs(pPosition[2])); + const float r(std::sqrt(x * x + y * y)); + + if (r < rInner) + rInner = r; + + if (absoluteZ < zMin) + zMin = absoluteZ; + } + + if (this->PassesQualityCuts(pTrack, trackParameters)) + { + const pandora::CartesianVector &momentumAtDca(trackParameters.m_momentumAtDca.Get()); + const float pX(momentumAtDca.GetX()), pY(momentumAtDca.GetY()), pZ(momentumAtDca.GetZ()); + const float pT(std::sqrt(pX * pX + pY * pY)); + + const float zCutForNonVertexTracks(m_trackerInnerR * std::fabs(pZ / pT) + m_settings.m_zCutForNonVertexTracks); + const bool passRzQualityCuts((zMin < zCutForNonVertexTracks) && (rInner < m_trackerInnerR + m_settings.m_maxBarrelTrackerInnerRDistance)); + + const bool isV0(this->IsV0(pTrack)); + const bool isDaughter(this->IsDaughter(pTrack)); + + // Decide whether track can be associated with a pandora cluster and used to form a charged PFO + //if ((d0 < m_settings.m_d0TrackCut) && (z0 < m_settings.m_z0TrackCut) && (rInner < m_trackerInnerR + m_settings.m_maxBarrelTrackerInnerRDistance)) + // FIXME! AD: (rInner < m_trackerInnerR + m_settings.m_maxBarrelTrackerInnerRDistance) check is removed + if ((d0 < m_settings.m_d0TrackCut) && (z0 < m_settings.m_z0TrackCut)) + { + canFormPfo = true; + } + else if (passRzQualityCuts && (0 != m_settings.m_usingNonVertexTracks)) + { + canFormPfo = true; + } + else if (isV0 || isDaughter) + { + canFormPfo = true; + } + + // Decide whether track can be used to form a charged PFO, even if track fails to be associated with a pandora cluster + const float particleMass(trackParameters.m_mass.Get()); + const float trackEnergy(std::sqrt(momentumAtDca.GetMagnitudeSquared() + particleMass * particleMass)); + + if ((0 != m_settings.m_usingUnmatchedVertexTracks) && (trackEnergy < m_settings.m_unmatchedVertexTrackMaxEnergy)) + { + if ((d0 < m_settings.m_d0UnmatchedVertexTrackCut) && (z0 < m_settings.m_z0UnmatchedVertexTrackCut) && + (rInner < m_trackerInnerR + m_settings.m_maxBarrelTrackerInnerRDistance)) + { + canFormClusterlessPfo = true; + } + else if (passRzQualityCuts && (0 != m_settings.m_usingNonVertexTracks) && (0 != m_settings.m_usingUnmatchedNonVertexTracks)) + { + canFormClusterlessPfo = true; + } + else if (isV0 || isDaughter) + { + canFormClusterlessPfo = true; + } + } + } + else if (this->IsDaughter(pTrack) || this->IsV0(pTrack)) + { + streamlog_out(WARNING) << "Recovering daughter or v0 track " << trackParameters.m_momentumAtDca.Get().GetMagnitude() << std::endl; + canFormPfo = true; + } + } + + trackParameters.m_canFormPfo = canFormPfo; + trackParameters.m_canFormClusterlessPfo = canFormClusterlessPfo; +} + +//------------------------------------------------------------------------------------------------------------------------------------------ + +void DDTrackCreatorALLEGRO::TrackReachesECAL(const EVENT::Track *const pTrack, PandoraApi::Track::Parameters &trackParameters) const +{ + // FIXME! AD: since currently we have do not have track hits -> check the radius of the trackState AtCalo. + // currently check is done only for ECAL barrel but should check if track reaches the Endcap! + pandora::CartesianVector posAtCalo(trackParameters.m_trackStateAtCalorimeter.Get().GetPosition()); + double radiusAtCalo = std::sqrt(posAtCalo.GetX()*posAtCalo.GetX() + posAtCalo.GetY()*posAtCalo.GetY()); + if(radiusAtCalo >= m_settings.m_eCalBarrelInnerR) trackParameters.m_reachesCalorimeter = true; + return; + + + // Calculate hit position information + float hitZMin(std::numeric_limits::max()); + float hitZMax(-std::numeric_limits::max()); + float hitOuterR(-std::numeric_limits::max()); + + int nTrackerHits(0); + int nEndcapDiskHits(0); + int maxOccupiedEndcapDiskLayer(0); + + const EVENT::TrackerHitVec &trackerHitVec(pTrack->getTrackerHits()); + const unsigned int nTrackHits(trackerHitVec.size()); + + for (unsigned int i = 0; i < nTrackHits; ++i) + { + const float x(static_cast(trackerHitVec[i]->getPosition()[0])); + const float y(static_cast(trackerHitVec[i]->getPosition()[1])); + const float z(static_cast(trackerHitVec[i]->getPosition()[2])); + const float r(std::sqrt(x * x + y * y)); + + if (z > hitZMax) + hitZMax = z; + + if (z < hitZMin) + hitZMin = z; + + if (r > hitOuterR) + hitOuterR = r; + + if ((r > m_trackerInnerR) && (r < m_trackerOuterR) && (std::fabs(z) <= m_trackerZmax)) + { + nTrackerHits++; + continue; + } + + for (unsigned int j = 0; j < m_nEndcapDiskLayers; ++j) + { + if ((r > m_endcapDiskInnerRadii[j]) && (r < m_endcapDiskOuterRadii[j]) && + (std::fabs(z) - m_settings.m_reachesECalFtdZMaxDistance < m_endcapDiskZPositions[j]) && + (std::fabs(z) + m_settings.m_reachesECalFtdZMaxDistance > m_endcapDiskZPositions[j])) + { + if (static_cast(j) > maxOccupiedEndcapDiskLayer) + maxOccupiedEndcapDiskLayer = static_cast(j); + + nEndcapDiskHits++; + break; + } + } + } + + // Require sufficient hits in barrel or endcap trackers, then compare extremal hit positions with tracker dimensions + if ((nTrackerHits >= m_settings.m_reachesECalNBarrelTrackerHits) || (nEndcapDiskHits >= m_settings.m_reachesECalNFtdHits)) + { + if ((hitOuterR - m_trackerOuterR > m_settings.m_reachesECalBarrelTrackerOuterDistance) || + (std::fabs(hitZMax) - m_trackerZmax > m_settings.m_reachesECalBarrelTrackerZMaxDistance) || + (std::fabs(hitZMin) - m_trackerZmax > m_settings.m_reachesECalBarrelTrackerZMaxDistance) || + (maxOccupiedEndcapDiskLayer >= m_settings.m_reachesECalMinFtdLayer)) + { + trackParameters.m_reachesCalorimeter = true; + return; + } + } + + // If track is lowpt, it may curl up and end inside tpc inner radius + const pandora::CartesianVector &momentumAtDca(trackParameters.m_momentumAtDca.Get()); + const float cosAngleAtDca(std::fabs(momentumAtDca.GetZ()) / momentumAtDca.GetMagnitude()); + const float pX(momentumAtDca.GetX()), pY(momentumAtDca.GetY()); + const float pT(std::sqrt(pX * pX + pY * pY)); + + if ((cosAngleAtDca > m_cosTracker) || (pT < m_settings.m_curvatureToMomentumFactor * m_settings.m_bField * m_trackerOuterR)) + { + trackParameters.m_reachesCalorimeter = true; + return; + } + + trackParameters.m_reachesCalorimeter = false; +} +